第二章 应力和应变
地震波传播的任何定量的描述,都要求其能表述固体介质的内力和变形的特征。现在我们对后面几章所需要的应力、应变理论的有关部分作简要的复习。虽然我们把这章作为独立的分析,但不对许多方程进行推导,读者想进一步了解其细节,可查阅连续介质力学的教科书。
三维介质的变形称为应变,介质不同部分之间的内力称为应力。应力和应变不是独立存在的,它们通过描述弹性固体性质的本构关系相联系。
2.1 应力的表述——应力张量
2.1.1应力表示
考虑一个在静力平衡状态下,均匀弹性介质里一个任意取向的无限小平面。
ˆ来规定。在nˆ方 平面的取向可以用这个平面的单位法向矢量n
向的一侧施加在此面单位面积上的力叫做牵引力,用矢量
ˆ相反方向的另一侧施加在此面上ˆ)=(tx,ty,tz)表示。在nt(n
ˆ)=-t(nˆ)。t在垂直的力与其大小相等,方向相反,即t(-n
于平面方向的分量叫做法应力,平行于平面方向的分量叫做
ˆ,这里P剪应力。在流体的情况下,没有剪应力,t=-pn
是压强。
上面的表示这是一个平面上的应力状况,为表示固体内部任意平面上的应力状态,应力张量τ在笛卡尔坐标系(图2.1)里可以用作用于yz,xz,xy平面的牵引力来定义+:
ˆ)tx(yˆ)tx(z⎡tx(xˆ)⎤⎡τxxτxyτxz⎤
⎥⎢⎥⎢
ˆ)ty(yˆ)ty(zˆ)⎥=⎢τyxτyyτyz⎥ (2.1) τ=⎢ty(x
⎢⎥⎢t(xˆ)tz(zˆ)⎥⎣zˆ)tz(y⎦⎣τzxτzyτzz⎦
在右式的表示中,第一个下角标表示面的法线方向,第二个下角标表示该面上应力在该坐标轴上的投影。
图2.1 在笛卡尔坐标系里描述作用在无限小立方体面上的力的牵引力矢
ˆ),t(yˆ),t(zˆ)。 量t(x
应力分量的符号规定如下:对于正应力,我们规定拉应力为正,压应力为负。对于剪应力,如果截面的外法线方向与坐标轴一致,则沿着坐标轴的正方向为正,反之为负;如果截面方向与外法线方向相反,则沿着坐标轴反方向为正。
下面讨论下角标颠倒后与颠倒前的值的关系。我们考虑xz面上在y轴方向延伸单位长度1的小微元立方体,在z方向的边长为∆z,在x方向的边长为∆x,如图所示,右边的τxz的外法线方向与x轴一致,因此沿z方向为正,而左边的τxz的外法线方向与x轴相反,逆z方向为正。τzx的分析与此分析类似。绕y轴的顺时针转动力矩为2τzxS
2τxzS
∆z∆z
=2τzx(1∆x)=τzx∆z∆x,逆时针旋转的力矩为22
∆x∆x
=2τxz(1∆z)=τxz∆z∆x,由于弹性体内部的微元不可能发生转动,因22
此两者必须相等,因此有τxz=τzx。类似地有:τxy=τyx,τyz=τzy。故应力张量是对称的,即:
⎡τxxτxyτxz⎤
⎢⎥
τ=τT=⎢τxyτyyτyz⎥ (2.2)
⎢τττ⎥⎣xzyzzz⎦
应力张量只包含6个独立的要素,它们足以完全描述介质中一个给定点的应力状态。
2、任意一个面上的应力可以由应力张量表示
一点的九个应力分量如果能够完全确定一点的应力状态,则其必须能够表达通过该点的任意斜截面上的应力矢量。
为了说明这一问题,在O点用三个坐标面和一任意斜截面截取一个微分四
面体单元,斜截面的法线方向矢量为n,它的三个方向余弦分别为l,m和n。
设斜截面上的应力为pn,i,j 和 k 分别为三个坐标轴方向的单位矢量,pn在坐标轴上的投影分别为px, py, pz。则应力矢量可以表示为
pn = pxi+ py j+ pz k
设S为ΔABC的面积,则
ΔOBC=lS, ΔOCA=mS, ΔOAB=nS
ΔABC的法线方向的单位矢量可表示为
n = l i+ m j + n k
微分四面体在应力矢量和体积力作用下应满足平衡条件。由x方向的平衡,可得
∑F
x
=0,
pxS-τxx∆OBC-τxy∆OAC-τxz∆OAB=0
注意,τxx,τxy,τxz取负是因为外法线方向与作用面的方向相反。将
公式
代入上式,有
pxS-τxxSl-τxySm-τxzSn=0
从而 px=τxxl+τxym+τxzn
py=τyxl+τyym+τyznpz=τzxl+τzym+τzzn
同理
如果采用张量记号,则上述公式可以表示为
上式给出了物体内一点的9个应力分量和通过同一点的各个微分面上的应力之间的关系。这一关系式表明,只要有了应力分量,就能够确定一点任意截面的应力矢量,或者正应力和切应力。因此应力分量可以确定一点的应力状态。
ˆ定义的平面,任一个取向由n一个侧面作用于另一个侧面的牵引力为应力张ˆ的乘积,即: 量与n
ˆ)⎤⎡τxxτxyτxz⎤⎡nˆx⎤⎡tx(n
⎥⎢⎥⎢⎥⎢
ˆ)=τnˆ=⎢ty(nˆ)⎥=⎢τxyτyyτyz⎥⎢nˆy⎥ (2.3) t(n
⎥⎢ˆ⎥⎢t(n⎥⎢
z⎦⎣zˆ)⎦⎣τxzτyzτzz⎦⎣n
ˆ的平面和xy,xz,yz平面所围限的四面体(柯西四面体)这可以通过对由垂直于n
面上的力求和作出说明。
ˆ的牵引力矢量t的线性算子,简单地说,应力张量是给出相对于法向矢量n
从这个意义上来说,应力张量与任何特定的坐标系无关。在地震学中,我们几乎总是把应力张量写成笛卡尔几何学里的一个3⨯3的矩阵。注意到对称的要求,应
力张量的独立参数由9个减少为6个,呈现为最一般形式的二阶张量(标量为零阶张量,矢量为一阶张量,等等)。
应力张量通常随在物质里的位置而变化,它是作用在固体里每一点的无限小的面上的力的度量。应力只给出了这些面由一边作用于另一边的力的度量,计量标准是单位面积上的力。然而,可能有其他力(如重力)。这些力称为体力,计量标准是每单位体积或单位质量上的力。
2.1.3坐标变换
如果一处的应力状态为σ11=1,σ22=-1,而其他元素为零,将其旋转45度后的应力状态为
0⎛
σ'=-1
0 ⎝
-10⎫
00⎪
⎪00⎪⎭
不同坐标系下的应力状态表示
一点的应力分量不仅随着变形体中点的位置在改变,而且即使在同一点,由于截面的法线方向不同,截面上的应力也不一样。假设已知在坐标系Oxyz中,弹性体中的某点的应力状态表示为:
⎡τxxτxyτxz⎤⎢⎥τ=⎢τxyτyyτyz⎥
⎢τττ⎥⎣xzyzzz⎦
则对于新的坐标系Ox’y’z’,这点的各应力分量是多少就是本小节讨论的问题。设新旧坐标系的方向余弦表如下:
的应力可以表示为:⎡τxxτxyτxz⎤⎡l1⎤⎢⎥
ˆ')=τxˆ'=⎢τxyτyyτyz⎥⎢m1⎥ t(x⎢⎥
⎢τττ⎥⎢n⎥⎣xzyzzz⎦⎣1⎦
ˆ'方向即得到在新坐标系中的各分量,将应力矢此应力有三个分量,将其再投影到x
ˆ'点乘即可得到: 量与x
⎡τxxτxyτxz⎤⎡l1⎤
⎢⎥
ˆ')=[l1m1n1]⎢τxyτyyτyz⎥⎢m1⎥τx'x'=i't(x⎢⎥ ⎢τττ⎥⎢n⎥
⎣xzyzzz⎦⎣1⎦
=l12τxx+m12τyy+n12τzz+2l1m1τxy+2m1n1τyz+2n1l1τxz
⎡τxxτxyτxz⎤⎡l2⎤⎢⎥
ˆ')=[l2m2n2]⎢τxyτyyτyz⎥⎢m2⎥τy'y'=j't(y⎢⎥ ⎢τττ⎥⎢n⎥
⎣xzyzzz⎦⎣2⎦
222=l2τxx+m2τyy+n2τzz+2l2m2τxy+2m2n2τyz+2n2l2τxz⎡τxxτxyτxz⎤⎡l3⎤
⎢⎥⎥ˆ')=[l3m3n3]⎢τxyτyyτyz⎥⎢τz'z'=k't(zm3⎢⎥ ⎢τττ⎥⎢n⎥⎣xzyzzz⎦⎣3⎦
22
=l32τxx+m3τyy+n3τzz+2l3m3τxy+2m3n3τyz+2n3l3τxz
⎡τxxτxyτxz⎤⎡l1⎤
⎢⎥
ˆ')=[l2m2n2]⎢τxyτyyτyz⎥⎢m1⎥τx'y'=j't(x⎢⎥ ⎢τττ⎥⎢n⎥
⎣xzyzzz⎦⎣1⎦
=l1l2τxx+m1m2τyy+n1n2τzz+(l1m2+l2m1)τxy+(m1n2+m2n1)τyz+(l2n1+n2l1)τxz⎡τxxτxyτxz⎤⎡l2⎤⎢⎥
ˆ')=[l3m3n3]⎢τxyτyyτyz⎥⎢m2⎥τy'z'=k't(y⎢⎥ ⎢τττ⎥⎢n⎥
⎣xzyzzz⎦⎣2⎦
=l2l3τxx+m2m3τyy+n2n3τzz+(l2m3+l3m2)τxy+(m2n3+m3n2)τyz+(l2n3+n2l3)τxz
⎡τxxτxyτxz⎤⎡l1⎤⎢⎥
ˆ')=[l3m3n3]⎢τxyτyyτyz⎥⎢m1⎥τx'z'=k't(x⎢⎥
⎢τττ⎥⎢n⎥⎣xzyzzz⎦⎣1⎦
=l1l3τxx+m1m3τyy+n1n3τzz+(l3m1+l1m3)τxy+(m3n1+m1n3)τyz+(l1n3+n1l3)τxz
…
所以
n1⎤⎡τxxτxyτxz⎤⎡l1l2l3⎤⎢⎥⎢⎥=NTτN n2⎥τττmmm⎢⎥xyyyyz123⎥⎢⎥⎢⎥nnn⎥n3⎥23⎦⎦⎣τxzτyzτzz⎦⎢⎣1
由上面的分析可知一个斜面上的正应力可以表示为:
⎡τxxτxyτxz⎤⎡nˆx⎤⎢⎥⎢⎥
ˆ)=nˆxnˆynˆz⎢τxyτyyτyz⎥⎢nˆy⎥ σ(n
⎢τττ⎥⎢n⎥ˆxzyzzzz⎣⎣⎦⎦⎡l1
τi'j'=⎢⎢l2
⎢⎣l3
m1
m2m3
[]
其剪应力为:
ˆ)]2-σ(nˆ)2 τ2=[t(n
2.1.4主应力和应力主轴
ˆ,使得在垂直于nˆ的面上,没有对任何应力张量,总是可以找到一个方向nˆ方向,在这种情况下: ˆ)沿n剪应力,也就是说,t(n
ˆ=λnˆ=τnˆt(n)
ˆ-λnˆ=0 (2.4) τnˆ=0(τ-Iλ)n
这里I是单位矩阵,λ是标量(不要把这些值同后面将讨论的拉梅参数相混淆)。这是一个本征值问题,只有
det[τ-Iλ]=0 (2.5)
τ-Iλ]=0的才有非零解。由于τ是对称的,是实数,所以本征值也是实数。det[
左端为λ的三次多项式,在线性代数中为矩阵的特征多项式。令其等于零,得到三个本征值λ1、λ2、λ3的解。将这三个解分别代入(2.5)就可以得到相应于三个
ˆ(1),nˆ(2),nˆ(3),它们是正交的,在弹性力学中称这本正值λ1、λ2、λ3的本征矢量为n
三个方向为应力主轴。垂直于应力主轴的平面叫做主平面。我们通过相似性变换,
ˆ(1),nˆ(2),nˆ(3) 的坐标系里: 把τ旋转到n
τR
⎡τ100⎤
⎥=NTτN=⎢0τ02⎢⎥
⎢⎣00τ3⎥⎦
(2.6)
这里τR是旋转的应力张量,τ1、τ2、τ3是主应力(与本征值λ1,λ2,λ3相等),N是本征矢量矩阵:
(1)⎡nx⎢N=⎢n(y1)
⎢(1)⎣nz
(2)nx
(3)
⎤nx⎥
n(y3)⎥ (2.7) (3)⎥nz⎦
n(y2)
(2)
nz
NT=N-1为归一化到单位长度的正交的本征矢量。
在MATLAB中矩阵的本征值和本征向量的求法为:
[X,D]=eig(sigma);其中sigma为三维应力矩阵,X为三个本征向量,D为本征值矩阵。得到本征向量和本征值后,你可以采用X’*sigma*X得到的结果验证是否是向量矩阵。
如果τ1=τ2=τ3,那么应力场处于流体静压状态,没有任何取向的面有剪应力。在流体的情况下,应力张量可写成:
⎡-P00⎤
⎥ (2.8) τ=⎢0-P0⎢⎥
⎢0-P⎥⎣0⎦
这里P是压强。
对于垂直应力不变的应力状态,其水平方向不同的应力状态可以表示为:
⎡τxxτxy0⎤⎡τxx-στxy0⎤⎢⎥⎢⎥
τ=⎢τxyτyy0⎥,其主应力方程可以表示为⎢τxyτyy-σ0⎥=0,或者展开为
⎢0⎢00τzz⎥0τzz-σ⎥⎣⎦⎣⎦
2(τzz-σ)[σ2-(τxx+τyy)σ+(τxxτyy-τxy)]=0,解这个方程可得到三个主应力为:
2⎣
1σ2=⎡τxx+τyy-
⎢⎣2
122⎤σ1=⎡τ+τ+τ-τ+4τxyxxyyxxyy⎢⎥
⎦
τ
2⎤-τ+4τxxyyxy⎥ ⎦
2
及σ3=τzz
2.1.1 应力值
应力以单位面积上的力为计量单位,在国际单位制(SI)中的单位是:
1帕斯卡(pa)=1牛顿米-2
回顾一下,1牛顿=1千克米秒-2=105达因。另一个普遍采用的力的单位是巴:
1巴(bar)=105pa 1千巴=108pa=100Mpa
1兆巴(Mbar)=1011pa=100Gpa(千兆帕)
如表2.1用参考模型PREM(Dziewonski和Anderson,1981)所给出的值所示,在地球里压力随深度快速增大。在400公里的深度,压力达到13.4Gpa,在核幔边界达136Gpa,在内核边界达329 Gpa。与此对比,月球中心的压力仅为4.8 Gpa,相当于在地球150公里深度所达到的值(Latham等人,1969),这是由于月球的质量小得多。
表2.1 地球内部压力与深度的关系
这些压力是地球内部的流体静压力。深部的剪切应力值小得多,包括与地幔对流有关的应力和由地震波传播所产生的动态应力。静态应力可能存在于地壳上部脆性部分。测量地壳剪应力是现代研究的课题。应力的量级是有争议的问题,地壳应力可能在100巴和1000巴(10-100 Mpa)之间,在接近活动断层的区域,应力有降低的趋势(活动断层的作用降低了应力)。
2.2 应变张量 2.2.1 位移场表示
现在让我们来考虑怎样描述连续介质里点的位置的变化。任何一点与参考时间t0的相对位置都可以用一个矢量场来描述,即位移场为:
u(r0)=r-r0 (2.9)
这里r是现在的位置,r0是参考点的位置。位移场是一个重要的概念,在这本书中经常涉及到。它是位置变化的绝对度量。与此相对照,应变是位移场相对变化
的局部度量,即位移场空间梯度的度量。应变与材料的变形或形状的变化有关,而与位置的绝对变化无关。例如,张应变是按照长度的相对变化来定义的。如果把一根100米长的细绳,一端固定,在另一端均匀地拉长到101米,那么沿细绳位移场从0变化到1米。依此,在绳的任何地方,应变场为0.01(1%)的常数。
考虑离开参考位置x0一个小的距离的某一点x的位移u=(ux,uy,uz),对u的每一个分量作泰勒级数展开,有:
∂u∂u∂u
ux=ux(x0)+xdx+xdy+xdz
∂x∂y∂z∂uy∂uy∂uy
uy=uy(x0)+dx+dy+dz
∂x∂y∂z∂u∂u∂u
uz=uz(x0)+zdx+zdy+zdz
∂x∂y∂z
写成矩阵形式得到:
⎡∂ux
⎢∂x⎢⎡ux⎤
⎢∂uy⎢⎥
u(x)=⎢uy⎥=u(x0)+⎢
⎢∂x⎢u⎥
⎣z⎦⎢∂uz
⎢⎣∂x
∂ux∂y∂uy∂y∂uz∂y
∂ux⎤
⎡dx⎤∂z⎥⎢⎥⎥
∂uy⎥⎢⎥
⎥⎢dy⎥=u(x0)+Jd (2.10) ∂z⎥⎢⎥∂uz⎥⎢⎥
dz⎥⎥⎢⎣⎦∂z⎦
这里d=x-x0,假定偏导数
∂ux∂uy
,等很小,它们的乘积可以忽略。那么,根 ∂x∂x
据无限小应变理论,我们就可以忽略展开式中的高阶项。庆幸的是在地震学中地球的应变几乎总是很小的,以致于这种近似是恰当的。我们可以通过把J 分成对称和反对称部分,把刚性旋转部分分离出来:
⎡∂ux∂ux∂ux⎤⎢∂x∂y∂z⎥⎢⎥⎢∂uy∂uy∂uy⎥J=⎢⎥=e+Ω (2.11)
∂x∂y∂z⎢⎥⎢∂uz∂uz∂uz⎥⎢⎥∂x∂y∂z⎣⎦这里应变张量e是对称的(eij=eji),可表达为:
⎡∂ux⎢
∂x⎢
⎢1∂uy∂uxe=⎢(+)
1∂ux∂uy
(+)2∂y∂x
∂uy
1∂ux∂uz⎤
(+)⎥2∂z∂x⎥1∂uy∂uz⎥(+)⎥ (2.12) ⎢2∂x
∂y∂y⎢⎢1(∂uz+∂ux1∂uz∂uy⎢⎣2∂x
∂z)
2(∂y+∂z
)这里旋转张量Ω是反对称的,可表达为:
⎡⎢
01∂ux∂uy
⎢
2(∂y-∂x)Ω=⎢⎢1∂u⎢
-x∂uy
2(∂y-∂x)0
⎢⎢-1(∂ux-∂uz1∂u
∂u⎢2∂z
)-(y-z)
⎣∂x2∂z∂y读者可对e+Ω=J作检验。
2.2.2物理解释
2∂z∂y⎥
∂uz
⎥∂z⎥
⎥⎦
1∂ux∂uz2(∂z-∂x)⎤⎥⎥12(∂uy∂z-∂uz∂y)⎥⎥ ⎥
⎥0
⎥⎥⎦
2.13) (
为简单起见,考虑平面的几何问题来讨论形变分量和位移之间的关系。如图,经过弹性体内的任意一点P,沿x轴和y轴的正方向取两个微小长度的线段PA=dx,PB=dy。假定弹性体受力后,P,A,B分别移到P’,A’,B’。首先考察线段PA
∂u
和PB的线应变,设P在x方向的位移为u,则A点在x方向的位移将是u+dx,
∂x在y方向的位移为v+
∂v∂udx;B点在x方向的位移将是u+dy,在y方向的位移∂x∂y
为v+
∂u∂u∂v
dy。则PA线段的线应变为其伸长增加量dx与其长度dx之比,为,
∂x∂x∂y
同理,PB线段的线应变为
∂v
。 ∂y
对于上述二维情况,若以PA,PB看做是矩形的两个边,其面积为dxdy,则变形后的P’A’
的长度为dx+
∂u∂v
dx(忽略其他高阶小量),P’B’的长度为dy+dy(忽略其他高阶小量),∂x∂y
则小体积元的面积趋于零时,面积增量与原来面积的极限为:
∂u⎫⎛∂v⎫⎛dx+dxdy+dy⎪-dxdy ⎪
∂x⎭ ∂yS-S0⎝∂u∂v⎝⎭ =≈+
S0dxdy∂x∂y
同样的,对于三维情况,可以得到体积增量和原来体积之比的极限为:
∂u⎫⎛⎛dx+d⎪x
∂x⎭ V-V0⎝⎝=
V0
∂v⎫⎛∂w⎫d+dy+⎪-dz ∂y⎪∂z⎭⎭⎝dxdydz
dxdydz
∂u∂v∂w≈++=divu ∂x∂y∂z
因此(2-12)式对角线元素之和为∆=
∂ux∂uy∂uz
++=tr[e]=∇⋅u为体应变。 ∂x∂y∂z
在线性代数中表示为tr[e]=e11+e22+e33是e的迹。注意,膨胀是由位移场的散度给出的。
现在考察线段PA与PB之间的直角的改变,该改变量有两部分贡献:一部分为由y方向的位移v引起的,即x方向的线段PA的转角α,另一部分为由x方向的位移u引起的,即y方向的线段PB的转角β.由图的关系不难看出:
α≈tgα=
∂v∂u,β≈tgβ= 。在三维情况下也具有同样道理。因此e矩阵中的对角线元∂x∂y
素为线应变,非对角线元素为线段偏转角度之和的一半(采用弧度表示)。
Fig 3
参考上图,一个质点P﹙y、z﹚向逆时针方向扭转到P’﹙y、z﹚,扭转角度为ωx,若其扭转半径为r,根据几何关系可得到:
y
y=rcosα,z=rsinα
其位移形变为
v=-rωxsinα=-zωx w=rωxcosα=yωx
将其分別对y及z微分且相加,得出
ωx=
1⎛∂w∂v⎫
-⎪
2⎝∂y∂z⎭
1⎛∂u∂w⎫1⎛∂v∂u⎫
--⎪,这正是(2-13)式的几个元素 ⎪和ωz= 2⎝∂z∂x⎭2⎝∂x∂y⎭
同理得到ωy=
另外,我们在高等数学中学过旋度的表示为
i
rotu=∇⨯u=
∂∂xux
j∂∂yuy
k∂∂zuz
⎛∂u∂uy⎫⎛∂ux∂uz⎫⎛∂uy∂ux⎫= z---⎪i+ ⎪k⎪j+
∂y∂z∂z∂x∂x∂y⎭⎝⎝⎭⎝⎭
⎛∂u∂uy⎫⎛∂ux∂uz
因此Ω表示了作为刚体的旋转,其中 z-为以x轴为轴的旋转量,-⎪ ∂z⎭∂x⎝∂z⎝∂y
⎫
⎪为⎭
⎛∂uy∂ux⎫
-以y轴为轴的旋转量, ⎪为以z轴为轴的旋转量。因此(2-13)就代表刚度的旋⎝∂x∂y⎭
转。
下面以二维应变为例列举几种应变状态,以此可以理解不同应变状态下的变形。
图
二维情况下不同应变状态的变形举例
2.3 本征值和本征向量
在上式可以看到,像应力张量那样,应变张量是对称的,包含6个独立的参
ˆ,可以找到应变主轴,即:
数。通过计算位移的取向n
ˆ=enˆ
(2.18) u=λn
这里的意义可以仿照应力张量的形式自然得知,在此不再赘述。
在下节我们将发现用脚标符号有助于表达应变张量。方程(2.12)可写成:
1
eij=(
∂iuj+∂jui) (2.20)
2这里假定i和j的顺序为1至3(按x,y,z的取向)。我们用了符号∂xuy=
2.2.1 应变值
∂uy∂x
。
由于应变表示长度的变化除以长度,故应变是无量纲的。与远场地震波通过相联系的动态应变,基本上小于10-6。
2.3 线性的应力—应变关系
在弹性介质中,应力和应变通过应力—应变的本构关系联系起来。应力张量和应变张量之间最一般的线性关系可以写成:
τij=cijklekl=
k=1,3l=1,3
∑∑c
ijkekl
e (2.21)
这里cijkl叫做弹性张量。在此,我们开始采用脚标符号求和的惯用做法。在乘积中任何重复的脚标都意味着求和是从脚标1到3。方程 (2.21)假定是完全弹性的,当应力作用,材料发生变形时,能量没有损失和衰减(有时,这些效应可由。在第六章之前,我们不考虑非弹性状态和衰减。 cijkl为复数来模拟)
应力张量是一个有81个(34)元素的四阶张量。
然而,由于应力张量和应变张量的对称性,它们各自只有六个分量,此时,Cijkl = Cjikl = Cijlk = Cjilk。因此如果不讨论预应力问题,应力分量和应变分量之间的表达式可以表示为:
⎧τ11=c1111e11+c1122e22+c1133e33+c1112e12+c1123e23+c1113e13⎪τ=ce+ce+ce+ce+ce+ce
[***********][***********]⎪22
⎪⎪τ33=c3311e11+c3322e22+c3333e33+c3312e12+c3323e23+c3313e13
(2.22) ⎨
τ=ce+ce+ce+ce+ce+ce[***********][***********]⎪12
⎪τ23=c2311e11+c2322e22+c2333e33+c2312e12+c2323e23+c2313e13⎪⎪⎩τ13=c1311e11+c1322e22+c1333e33+c1312e12+c1323e23+c1313e13
上式只有36个弹性常数。可以证明,对于一个保守系统(即无能量损失),
C1122=C2211,C1133=C3311,C2233=C3322,C1213=C1312,C1223=C2312,C1323=C2313,C1121=C2111,C1113=C1311,C1123=C2311,C2212=C1222,C2213=C1322,C2223=C2322
,
因此对于极端的各向异性介质,C3321=C2133,C3313=C1333,C3323=C2333 (证明从略)。
独立的弹性参数为21个,这些元素只有21个是独立的。这21个元素是确定弹
性固体的最一般形式的应力—应变关系所必须的。
这样一种固体的性质可能随方向变化,如果是这样的话,就称这种介质是各向异性的。与此相反,各向同性的固体,所有方向的性质是相同的。对地球内部,大多数情况下,各向同性是合适的一级近似。但在一些地区观测到各向异性,这是现代研究的一个重要领域。
如果我们作了各向同性的假定,独立参数减至2个:λ和μ,这些参数又叫做介质的拉梅参数。
Cijkl=λδijδkl+μ(δikδjl+δilδjk)
可
见
,
C1
1
=
1
C1=C=μ2λ222
,+
C1122=C1133=C2233=λ
,,,
C1313=C1212=C2323=μ
C1112=C1113=C1123=C2212=C2213=C2223=C3312=C3313=C3323=0C1213=C1223=C1323=0。
我们将看到,拉梅参数及介质的密度将最终确定介质的地震波速度。对各向同性固体,应力—应变方程(2.21)为:
τij=⎡⎣λδijδkl+μ(δilδjk+δikδjl)⎤⎦ekl=λδijekk+2μeij (2.23)
这里我们用eij=eji合并μ的项。注意ekk=tr[e]是e的对角线元素的和。用这方程,我们根据应变可直接写出应力张量:
⎡λtr[e]+2μe112μe12⎢
τ=⎢2μe21λtr[e]+2μe22
⎢2μe2μe32⎣31
⎤⎥
2μe23⎥ (2.24) λtr[e]+2μe33⎥⎦2μe13
μ叫做剪这两个拉梅参数完全描述了各向同性固体里的线性的应力—应变关系。切模量,是介质抗剪切的度量。它的值是所作用的剪切应力与所导致的剪应变比率的一半,即μ=τxy/2exy(这里的2与应变张量的1/2有关)。另一个拉梅参数
λ,没有简单的物理解释。各向同性固体的其他普遍使用的弹性参数中很多具有
特定的物理意义。
2.4 弹性参数的物理意义及相互转换
杨氏模量E:两端拉伸的柱体的张应力与张应变的比率,
正应力和线应变之间满足胡克定律:
σε=
E
所以E表示弹性材料抵抗拉张(或压缩)的能力,E越大说明弹性体越难变形。 Possion比
弹性体纵向伸长(或缩短)∆L后,其横向也会有∆d的变化,
∆d
γ=-
叫做Possion系数。式中负号表示纵向应变和横向应变方向相反,为保证γ为正而取的。Possion系数仅与材料本身的性质有关。实验表明,对于一切介质,γ介于0到0.5之间,对于地幔介质,常用0.25表示地幔的大部分,对于地球外核(液态),取为0.5.
体积模量K:在实际地球中,只受单向压力或拉力的情形很少,一般情况下是各个方向都受力,最常见的是液体静压力。弹性体在静压力P作用下提及变化∆V,则有关系
∆V
P=-K
V
因此体变模量K为流体静压力与由其所导致的体积变化的比率,是介质不可压缩性的度量。
应指出,体变模量K可以从杨氏模量E和Possion系数推导出来。设一个立方体边长为L,体积则为V=L3,在液体静压应力下,整个体积的相对变化率可由
∆V3∆L
=边长的相对变化率决定,即:,就任何一个方向而言,在这个方向上VL
的压应力,使其长度改变∆L'而在与此方向垂直的另两个方向的压应力,又使其长度改变了-2γ∆L',故总的效果是这三个方向压应力作用之和:∆L=(1-2γ)∆L',从而得出:
∆V3∆L∆L'
==3(1-2γ),注意,∆L'为仅在单向压力作用下的长度的改VLL
变量,满足关系式:
∆L'P∆L'=,P=E LEL这样可得到:
ε=
K=
E
31-2γ弹性系数之间的关系: 单向拉伸实验:
τ11=σ,τ22=τ33=τ23=τ13=τ12=0
E
由胡克定律可得: e11=
τ11
,e22=e33=-γe11=-γ
τ11
E
τ11=2μe11+λekk0=2μe22+λekk 0=2μe33+λekk
三式相加得:
ekk=τ11 3λ+2μμ(3λ+2μ) λ+μ代入上式第一式得:E=
代入第二式或第三式得:γ=
各向均匀压缩试验: λ 2(λ+μ)
τ11=τ22=τ33=-p,τ12=τ13=τ23=0 体变模量定义为:K=
由广义胡克定律: -p ekk
-p=2μe11+λekk
-p=2μe22+λekk
-p=2μe33+λekk
三式相加得:
-3p=(3λ+2μ)ekk 与体变模量的定义比较可得K=λ+
在地震学中,我们主要涉及压缩(P)和剪切(S)波的速度。如后面第三章所指出,它们可由这些弹性常数和密度ρ来计算:
P波的速度α表达为: 2μ (2.26) 3
α=
S波的速度β表达为: +2 (2.28) ρ
β= (2.29) ρ
泊松比往往用来度量P波和S波速度的相对大小,表达为:
α2-2β2
(2.30) γ=222(α-β)
注意γ无量纲,在0-0.5之间变化,下限相当于流体的情况(μ=0)。对泊松固体,λ=μ,σ=0.25,α/β=3。多数地壳岩石,泊松比在0.25与0.30之间。
2.3.1 弹性模量的单位
拉梅参数、扬氏模量、体积模量都有像应力那样的单位(帕斯卡)。回顾一下:
1帕斯卡(pa)=1牛顿米-2=1千克米-1秒-2
注意当其被密度除时,结果即为速度平方(适合于方程2.28和2.29)的单位。
练习:
0-1⎫⎛1 ⎪1、地壳中某点的应力状态在北东下坐标系中表示为:τij= 0-10⎪,
-101⎪⎝⎭
单位MPa。求其三个主应力及其方向。
2、已知地壳中某点的主应力为σ1=75bar, σ2=50bar, σ3=-50bar,一斜截面的法线与三个主轴成等角,求该面上的正应力和剪应力。(提示该面的方向余弦为1
,1
,1
3)。
⎛111⎫ ⎪3、已知某点的应力张量为:τ= 111⎪,求其主应力大小和主轴方向。
111⎪⎝⎭
2.1 用方程(2.4),(2.18)和(2.24)说明对各向同性介质,应力主轴总是恰好与应变主轴重合。
2.2 由方程(2.28)和(2.29),根据地震波速度和密度导出拉梅参数的表达式。
2.3对于均匀各向同性弹性体,用应变分量表示应力分量的胡克定律τij=λδijekk+2μeij,试利用弹性常数之间的关系导出应力分量表示为应变分量的表达式:eij=-γ
Eδijτkk+1+γτij。 E
2.3 对泊松比为0.30的岩石,P/S波的速度比是多少?
2.4 观测到实验室里的花岗岩样品的P波速度为5.5km/s,密度为2.6Mg/m3。假定样品是泊松固体,求拉梅参数,扬氏模量、体积模量的值。给出你的以帕斯
卡为单位的答案。如果样品以地球24公里深度所存在的压力压缩,那么这样品体积的相对变化是多少?对这问题,假定当样品压缩时,体积模量没有变化。
2.5 用PREM模型的值(附录1)计算(a)核幔边界(CMB)和(b)内核边界(ICB)两边的体积模量的值。以帕斯卡单位表达你的答案。
2.6 圣地亚哥,加利福尼亚大学在圣地亚哥东北山脉(在Anza附近)建立了Pinon Flat观测台(PFO),仪器包括测量地壳变形的高质量的应变计。
(a) 假定在PFO底下5公里,地震波的速度α=6km/s,β=3.5km/s,密度ρ=2.7Mg/m3,根据这些参数计算拉梅参数λ,μ的值。以帕斯卡单位给出答案。
(b)1922年南加利福尼亚PFO以北80公里的Landers地震(Ms=7.3)后,PFO应变仪观测到相对于地震前,应变有一个大的静态变化。应变张量的水平分量按以下数量变化:e11=-0.26⨯10-6,e22=0.92⨯10-6,e12=-0.69⨯10-6,这里脚标1表示东,2表示北,拉张为正。你可以假定应变的变化是在地震时瞬间出现的。假定在深部,这些值也是正确的,用你在(a)所得到的结果来确定由于地震,在5公里深的地方造成的应力变化,即计算τ11,τ22和τ12的变化。把这看成是假定在垂直方向应变没有变化,应变与深度无关的二维问题。
(c)计算在PFO因Lander地震,应变主轴的取向(水平)。用方位角表达你的答案(北东多少度)。
(d)在PFO观测到应变稳定的长期变化,一年的变化量为:
注意到长期的应变变化接近e11=0.101⨯10-6,e22=-0.02⨯10-6,e12=0.005⨯10-6。
于简单的E-W向拉伸,假定在过去100年,应变的速率是稳定的,初始应力为零,计算5公里深度应力张量的分量(注意,这或许是不很真实的假设!)不包括在5公里深度应力的流体静压分量。
图2.4 在1992年南加利福尼亚地震(MS 7.3)震中
以南80公里的PFO观测台观测到的应变变化。
(e)农民Bob在PFO附近有一平方公里的一小片土地,他将其围起来,并作高精度的测量。农民Bob每年增加或损失多少土地?由于Landers地震,他增加或损失多少土地?
(f)(计算)编写计算机程序计算垂直断层两侧从0°到170°之间不同方位(从北到东以10°增加)的应力。对你在(b)和(d)中计算的应力张量,作表列出断层方位和相应的剪应力及断层两边的法应力(对Landers地震,这些值是应力的变化,不是绝对值)。对每种情况,什么方位有最大的剪应力?
对于(b)中的应力状态:
Vp=6000;Vs=3500;dens=2700;
miu=dens*Vs^2;
lam=dens*Vp^2-2*miu;
e11=0.26e-6;e22=0.92e-6;e12=-0.69e-6;
ekk=e11+e22;
s11=lam*(ekk)+2*miu*e11;
s12=2*miu*e12;
s22=lam*(ekk)+2*miu*e22;
s=[s11,s12;s12,s22];
for theta=0:10:170
N=[cos(deg2rad(theta)) -sin(deg2rad(theta));sin(deg2rad(theta))
cos((deg2rad(theta)))];
theta
S=N*s*N'
end
对于(d)的情况,改变应变值即可。
(g)(计算)一些研究(例如,Stein等人,1992,1994;Harris和Simpson,1992;Harris等人,1995)通过假定沿断层地震破裂的可能性与库伦破裂函数CFF有关,模拟了大地震的空间分布。忽略了孔隙流体压力效应,CFF的变化可以表达为:
∆CFF=∆s+μ∆τn
这里τs是剪应力(曳引力),τn是法应力,μ是静摩擦系数(不能同剪切模量混淆)。注意当剪应力增大,断层上的压应力减小时,CFF增大(回顾一下,在我们的符号约定中,张应力为正,压应力为负)。假定μ=0.2,修改你的计算机程序来计算每个断层取向的∆CFF。作出由于每个断层方位长期的变化,∆CFF年变化的表。
Vp=6000;Vs=3500;dens=2700;
miu=dens*Vs^2;
lam=dens*Vp^2-2*miu;
e11=0.26e-6;e22=0.92e-6;e12=-0.69e-6;
ekk=e11+e22;
s11=lam*(ekk)+2*miu*e11;
s12=2*miu*e12;
s22=lam*(ekk)+2*miu*e22;
s=[s11,s12;s12,s22];
for theta=0:10:170
N=[sin(deg2rad(theta)) cos(deg2rad(theta))];
t=s*N';
sigmaN=N*t;
shear=sqrt(sum(t.^2)-sigmaN^2);
theta
CFF=abs(shear)+0.2*sigmaN
end
(h)(计算)现在假定当他们的长期库伦破裂函数CFF达到某临界阀值时,断层将破裂。至下次地震的时间变量可表达为:
∆t=CFF1000+L-CFF1000 CFFa
这里CFFa是CFF的周年变化,CFF1000是CFF的千年变化,CFF1000+L是CFF的千年变化+Landers地震的变化(注意CFF1000+L≠CFF1000+CFFL)。用在下次地震前每年断层取向加速或减速的时间计算Landers地震的效应。用在下次地震前,加速为正的时间,减速为负的时间的符号约定,陈述你的每年的答案(敬告:这是复杂的)。把断层面上剪应力值的符号反过来,审核你的答案。通常(但不总是)当长期的和Landers地震的剪切变化符号一致(两者都是正或都是负)时,地震时间将提前,当剪切变化的符号不一致时,地震时间将推迟。
(i)提示:对(f)-(h),尤其是(h)作符号校正可能是复杂的。应力可能是正或负。为了得到正确的结果,对每个断层方位规定两个单位矢量,一个平
ˆ)ˆ计算牵引力矢量。行断层(fˆ),另一个垂直于断层(p。通过把应力张量乘以p
ˆ的点积把牵引力矢量分解为剪应力和法应力。然后通过分别计算fˆ与p自然,此
ˆ应是单位长度。 时fˆ和p
(j)Landers地震后,没有观测到PFO附近地震活动性(小震活动)增强。这说明临界值的CFF模型是正确的吗?
第二章 应力和应变
地震波传播的任何定量的描述,都要求其能表述固体介质的内力和变形的特征。现在我们对后面几章所需要的应力、应变理论的有关部分作简要的复习。虽然我们把这章作为独立的分析,但不对许多方程进行推导,读者想进一步了解其细节,可查阅连续介质力学的教科书。
三维介质的变形称为应变,介质不同部分之间的内力称为应力。应力和应变不是独立存在的,它们通过描述弹性固体性质的本构关系相联系。
2.1 应力的表述——应力张量
2.1.1应力表示
考虑一个在静力平衡状态下,均匀弹性介质里一个任意取向的无限小平面。
ˆ来规定。在nˆ方 平面的取向可以用这个平面的单位法向矢量n
向的一侧施加在此面单位面积上的力叫做牵引力,用矢量
ˆ相反方向的另一侧施加在此面上ˆ)=(tx,ty,tz)表示。在nt(n
ˆ)=-t(nˆ)。t在垂直的力与其大小相等,方向相反,即t(-n
于平面方向的分量叫做法应力,平行于平面方向的分量叫做
ˆ,这里P剪应力。在流体的情况下,没有剪应力,t=-pn
是压强。
上面的表示这是一个平面上的应力状况,为表示固体内部任意平面上的应力状态,应力张量τ在笛卡尔坐标系(图2.1)里可以用作用于yz,xz,xy平面的牵引力来定义+:
ˆ)tx(yˆ)tx(z⎡tx(xˆ)⎤⎡τxxτxyτxz⎤
⎥⎢⎥⎢
ˆ)ty(yˆ)ty(zˆ)⎥=⎢τyxτyyτyz⎥ (2.1) τ=⎢ty(x
⎢⎥⎢t(xˆ)tz(zˆ)⎥⎣zˆ)tz(y⎦⎣τzxτzyτzz⎦
在右式的表示中,第一个下角标表示面的法线方向,第二个下角标表示该面上应力在该坐标轴上的投影。
图2.1 在笛卡尔坐标系里描述作用在无限小立方体面上的力的牵引力矢
ˆ),t(yˆ),t(zˆ)。 量t(x
应力分量的符号规定如下:对于正应力,我们规定拉应力为正,压应力为负。对于剪应力,如果截面的外法线方向与坐标轴一致,则沿着坐标轴的正方向为正,反之为负;如果截面方向与外法线方向相反,则沿着坐标轴反方向为正。
下面讨论下角标颠倒后与颠倒前的值的关系。我们考虑xz面上在y轴方向延伸单位长度1的小微元立方体,在z方向的边长为∆z,在x方向的边长为∆x,如图所示,右边的τxz的外法线方向与x轴一致,因此沿z方向为正,而左边的τxz的外法线方向与x轴相反,逆z方向为正。τzx的分析与此分析类似。绕y轴的顺时针转动力矩为2τzxS
2τxzS
∆z∆z
=2τzx(1∆x)=τzx∆z∆x,逆时针旋转的力矩为22
∆x∆x
=2τxz(1∆z)=τxz∆z∆x,由于弹性体内部的微元不可能发生转动,因22
此两者必须相等,因此有τxz=τzx。类似地有:τxy=τyx,τyz=τzy。故应力张量是对称的,即:
⎡τxxτxyτxz⎤
⎢⎥
τ=τT=⎢τxyτyyτyz⎥ (2.2)
⎢τττ⎥⎣xzyzzz⎦
应力张量只包含6个独立的要素,它们足以完全描述介质中一个给定点的应力状态。
2、任意一个面上的应力可以由应力张量表示
一点的九个应力分量如果能够完全确定一点的应力状态,则其必须能够表达通过该点的任意斜截面上的应力矢量。
为了说明这一问题,在O点用三个坐标面和一任意斜截面截取一个微分四
面体单元,斜截面的法线方向矢量为n,它的三个方向余弦分别为l,m和n。
设斜截面上的应力为pn,i,j 和 k 分别为三个坐标轴方向的单位矢量,pn在坐标轴上的投影分别为px, py, pz。则应力矢量可以表示为
pn = pxi+ py j+ pz k
设S为ΔABC的面积,则
ΔOBC=lS, ΔOCA=mS, ΔOAB=nS
ΔABC的法线方向的单位矢量可表示为
n = l i+ m j + n k
微分四面体在应力矢量和体积力作用下应满足平衡条件。由x方向的平衡,可得
∑F
x
=0,
pxS-τxx∆OBC-τxy∆OAC-τxz∆OAB=0
注意,τxx,τxy,τxz取负是因为外法线方向与作用面的方向相反。将
公式
代入上式,有
pxS-τxxSl-τxySm-τxzSn=0
从而 px=τxxl+τxym+τxzn
py=τyxl+τyym+τyznpz=τzxl+τzym+τzzn
同理
如果采用张量记号,则上述公式可以表示为
上式给出了物体内一点的9个应力分量和通过同一点的各个微分面上的应力之间的关系。这一关系式表明,只要有了应力分量,就能够确定一点任意截面的应力矢量,或者正应力和切应力。因此应力分量可以确定一点的应力状态。
ˆ定义的平面,任一个取向由n一个侧面作用于另一个侧面的牵引力为应力张ˆ的乘积,即: 量与n
ˆ)⎤⎡τxxτxyτxz⎤⎡nˆx⎤⎡tx(n
⎥⎢⎥⎢⎥⎢
ˆ)=τnˆ=⎢ty(nˆ)⎥=⎢τxyτyyτyz⎥⎢nˆy⎥ (2.3) t(n
⎥⎢ˆ⎥⎢t(n⎥⎢
z⎦⎣zˆ)⎦⎣τxzτyzτzz⎦⎣n
ˆ的平面和xy,xz,yz平面所围限的四面体(柯西四面体)这可以通过对由垂直于n
面上的力求和作出说明。
ˆ的牵引力矢量t的线性算子,简单地说,应力张量是给出相对于法向矢量n
从这个意义上来说,应力张量与任何特定的坐标系无关。在地震学中,我们几乎总是把应力张量写成笛卡尔几何学里的一个3⨯3的矩阵。注意到对称的要求,应
力张量的独立参数由9个减少为6个,呈现为最一般形式的二阶张量(标量为零阶张量,矢量为一阶张量,等等)。
应力张量通常随在物质里的位置而变化,它是作用在固体里每一点的无限小的面上的力的度量。应力只给出了这些面由一边作用于另一边的力的度量,计量标准是单位面积上的力。然而,可能有其他力(如重力)。这些力称为体力,计量标准是每单位体积或单位质量上的力。
2.1.3坐标变换
如果一处的应力状态为σ11=1,σ22=-1,而其他元素为零,将其旋转45度后的应力状态为
0⎛
σ'=-1
0 ⎝
-10⎫
00⎪
⎪00⎪⎭
不同坐标系下的应力状态表示
一点的应力分量不仅随着变形体中点的位置在改变,而且即使在同一点,由于截面的法线方向不同,截面上的应力也不一样。假设已知在坐标系Oxyz中,弹性体中的某点的应力状态表示为:
⎡τxxτxyτxz⎤⎢⎥τ=⎢τxyτyyτyz⎥
⎢τττ⎥⎣xzyzzz⎦
则对于新的坐标系Ox’y’z’,这点的各应力分量是多少就是本小节讨论的问题。设新旧坐标系的方向余弦表如下:
的应力可以表示为:⎡τxxτxyτxz⎤⎡l1⎤⎢⎥
ˆ')=τxˆ'=⎢τxyτyyτyz⎥⎢m1⎥ t(x⎢⎥
⎢τττ⎥⎢n⎥⎣xzyzzz⎦⎣1⎦
ˆ'方向即得到在新坐标系中的各分量,将应力矢此应力有三个分量,将其再投影到x
ˆ'点乘即可得到: 量与x
⎡τxxτxyτxz⎤⎡l1⎤
⎢⎥
ˆ')=[l1m1n1]⎢τxyτyyτyz⎥⎢m1⎥τx'x'=i't(x⎢⎥ ⎢τττ⎥⎢n⎥
⎣xzyzzz⎦⎣1⎦
=l12τxx+m12τyy+n12τzz+2l1m1τxy+2m1n1τyz+2n1l1τxz
⎡τxxτxyτxz⎤⎡l2⎤⎢⎥
ˆ')=[l2m2n2]⎢τxyτyyτyz⎥⎢m2⎥τy'y'=j't(y⎢⎥ ⎢τττ⎥⎢n⎥
⎣xzyzzz⎦⎣2⎦
222=l2τxx+m2τyy+n2τzz+2l2m2τxy+2m2n2τyz+2n2l2τxz⎡τxxτxyτxz⎤⎡l3⎤
⎢⎥⎥ˆ')=[l3m3n3]⎢τxyτyyτyz⎥⎢τz'z'=k't(zm3⎢⎥ ⎢τττ⎥⎢n⎥⎣xzyzzz⎦⎣3⎦
22
=l32τxx+m3τyy+n3τzz+2l3m3τxy+2m3n3τyz+2n3l3τxz
⎡τxxτxyτxz⎤⎡l1⎤
⎢⎥
ˆ')=[l2m2n2]⎢τxyτyyτyz⎥⎢m1⎥τx'y'=j't(x⎢⎥ ⎢τττ⎥⎢n⎥
⎣xzyzzz⎦⎣1⎦
=l1l2τxx+m1m2τyy+n1n2τzz+(l1m2+l2m1)τxy+(m1n2+m2n1)τyz+(l2n1+n2l1)τxz⎡τxxτxyτxz⎤⎡l2⎤⎢⎥
ˆ')=[l3m3n3]⎢τxyτyyτyz⎥⎢m2⎥τy'z'=k't(y⎢⎥ ⎢τττ⎥⎢n⎥
⎣xzyzzz⎦⎣2⎦
=l2l3τxx+m2m3τyy+n2n3τzz+(l2m3+l3m2)τxy+(m2n3+m3n2)τyz+(l2n3+n2l3)τxz
⎡τxxτxyτxz⎤⎡l1⎤⎢⎥
ˆ')=[l3m3n3]⎢τxyτyyτyz⎥⎢m1⎥τx'z'=k't(x⎢⎥
⎢τττ⎥⎢n⎥⎣xzyzzz⎦⎣1⎦
=l1l3τxx+m1m3τyy+n1n3τzz+(l3m1+l1m3)τxy+(m3n1+m1n3)τyz+(l1n3+n1l3)τxz
…
所以
n1⎤⎡τxxτxyτxz⎤⎡l1l2l3⎤⎢⎥⎢⎥=NTτN n2⎥τττmmm⎢⎥xyyyyz123⎥⎢⎥⎢⎥nnn⎥n3⎥23⎦⎦⎣τxzτyzτzz⎦⎢⎣1
由上面的分析可知一个斜面上的正应力可以表示为:
⎡τxxτxyτxz⎤⎡nˆx⎤⎢⎥⎢⎥
ˆ)=nˆxnˆynˆz⎢τxyτyyτyz⎥⎢nˆy⎥ σ(n
⎢τττ⎥⎢n⎥ˆxzyzzzz⎣⎣⎦⎦⎡l1
τi'j'=⎢⎢l2
⎢⎣l3
m1
m2m3
[]
其剪应力为:
ˆ)]2-σ(nˆ)2 τ2=[t(n
2.1.4主应力和应力主轴
ˆ,使得在垂直于nˆ的面上,没有对任何应力张量,总是可以找到一个方向nˆ方向,在这种情况下: ˆ)沿n剪应力,也就是说,t(n
ˆ=λnˆ=τnˆt(n)
ˆ-λnˆ=0 (2.4) τnˆ=0(τ-Iλ)n
这里I是单位矩阵,λ是标量(不要把这些值同后面将讨论的拉梅参数相混淆)。这是一个本征值问题,只有
det[τ-Iλ]=0 (2.5)
τ-Iλ]=0的才有非零解。由于τ是对称的,是实数,所以本征值也是实数。det[
左端为λ的三次多项式,在线性代数中为矩阵的特征多项式。令其等于零,得到三个本征值λ1、λ2、λ3的解。将这三个解分别代入(2.5)就可以得到相应于三个
ˆ(1),nˆ(2),nˆ(3),它们是正交的,在弹性力学中称这本正值λ1、λ2、λ3的本征矢量为n
三个方向为应力主轴。垂直于应力主轴的平面叫做主平面。我们通过相似性变换,
ˆ(1),nˆ(2),nˆ(3) 的坐标系里: 把τ旋转到n
τR
⎡τ100⎤
⎥=NTτN=⎢0τ02⎢⎥
⎢⎣00τ3⎥⎦
(2.6)
这里τR是旋转的应力张量,τ1、τ2、τ3是主应力(与本征值λ1,λ2,λ3相等),N是本征矢量矩阵:
(1)⎡nx⎢N=⎢n(y1)
⎢(1)⎣nz
(2)nx
(3)
⎤nx⎥
n(y3)⎥ (2.7) (3)⎥nz⎦
n(y2)
(2)
nz
NT=N-1为归一化到单位长度的正交的本征矢量。
在MATLAB中矩阵的本征值和本征向量的求法为:
[X,D]=eig(sigma);其中sigma为三维应力矩阵,X为三个本征向量,D为本征值矩阵。得到本征向量和本征值后,你可以采用X’*sigma*X得到的结果验证是否是向量矩阵。
如果τ1=τ2=τ3,那么应力场处于流体静压状态,没有任何取向的面有剪应力。在流体的情况下,应力张量可写成:
⎡-P00⎤
⎥ (2.8) τ=⎢0-P0⎢⎥
⎢0-P⎥⎣0⎦
这里P是压强。
对于垂直应力不变的应力状态,其水平方向不同的应力状态可以表示为:
⎡τxxτxy0⎤⎡τxx-στxy0⎤⎢⎥⎢⎥
τ=⎢τxyτyy0⎥,其主应力方程可以表示为⎢τxyτyy-σ0⎥=0,或者展开为
⎢0⎢00τzz⎥0τzz-σ⎥⎣⎦⎣⎦
2(τzz-σ)[σ2-(τxx+τyy)σ+(τxxτyy-τxy)]=0,解这个方程可得到三个主应力为:
2⎣
1σ2=⎡τxx+τyy-
⎢⎣2
122⎤σ1=⎡τ+τ+τ-τ+4τxyxxyyxxyy⎢⎥
⎦
τ
2⎤-τ+4τxxyyxy⎥ ⎦
2
及σ3=τzz
2.1.1 应力值
应力以单位面积上的力为计量单位,在国际单位制(SI)中的单位是:
1帕斯卡(pa)=1牛顿米-2
回顾一下,1牛顿=1千克米秒-2=105达因。另一个普遍采用的力的单位是巴:
1巴(bar)=105pa 1千巴=108pa=100Mpa
1兆巴(Mbar)=1011pa=100Gpa(千兆帕)
如表2.1用参考模型PREM(Dziewonski和Anderson,1981)所给出的值所示,在地球里压力随深度快速增大。在400公里的深度,压力达到13.4Gpa,在核幔边界达136Gpa,在内核边界达329 Gpa。与此对比,月球中心的压力仅为4.8 Gpa,相当于在地球150公里深度所达到的值(Latham等人,1969),这是由于月球的质量小得多。
表2.1 地球内部压力与深度的关系
这些压力是地球内部的流体静压力。深部的剪切应力值小得多,包括与地幔对流有关的应力和由地震波传播所产生的动态应力。静态应力可能存在于地壳上部脆性部分。测量地壳剪应力是现代研究的课题。应力的量级是有争议的问题,地壳应力可能在100巴和1000巴(10-100 Mpa)之间,在接近活动断层的区域,应力有降低的趋势(活动断层的作用降低了应力)。
2.2 应变张量 2.2.1 位移场表示
现在让我们来考虑怎样描述连续介质里点的位置的变化。任何一点与参考时间t0的相对位置都可以用一个矢量场来描述,即位移场为:
u(r0)=r-r0 (2.9)
这里r是现在的位置,r0是参考点的位置。位移场是一个重要的概念,在这本书中经常涉及到。它是位置变化的绝对度量。与此相对照,应变是位移场相对变化
的局部度量,即位移场空间梯度的度量。应变与材料的变形或形状的变化有关,而与位置的绝对变化无关。例如,张应变是按照长度的相对变化来定义的。如果把一根100米长的细绳,一端固定,在另一端均匀地拉长到101米,那么沿细绳位移场从0变化到1米。依此,在绳的任何地方,应变场为0.01(1%)的常数。
考虑离开参考位置x0一个小的距离的某一点x的位移u=(ux,uy,uz),对u的每一个分量作泰勒级数展开,有:
∂u∂u∂u
ux=ux(x0)+xdx+xdy+xdz
∂x∂y∂z∂uy∂uy∂uy
uy=uy(x0)+dx+dy+dz
∂x∂y∂z∂u∂u∂u
uz=uz(x0)+zdx+zdy+zdz
∂x∂y∂z
写成矩阵形式得到:
⎡∂ux
⎢∂x⎢⎡ux⎤
⎢∂uy⎢⎥
u(x)=⎢uy⎥=u(x0)+⎢
⎢∂x⎢u⎥
⎣z⎦⎢∂uz
⎢⎣∂x
∂ux∂y∂uy∂y∂uz∂y
∂ux⎤
⎡dx⎤∂z⎥⎢⎥⎥
∂uy⎥⎢⎥
⎥⎢dy⎥=u(x0)+Jd (2.10) ∂z⎥⎢⎥∂uz⎥⎢⎥
dz⎥⎥⎢⎣⎦∂z⎦
这里d=x-x0,假定偏导数
∂ux∂uy
,等很小,它们的乘积可以忽略。那么,根 ∂x∂x
据无限小应变理论,我们就可以忽略展开式中的高阶项。庆幸的是在地震学中地球的应变几乎总是很小的,以致于这种近似是恰当的。我们可以通过把J 分成对称和反对称部分,把刚性旋转部分分离出来:
⎡∂ux∂ux∂ux⎤⎢∂x∂y∂z⎥⎢⎥⎢∂uy∂uy∂uy⎥J=⎢⎥=e+Ω (2.11)
∂x∂y∂z⎢⎥⎢∂uz∂uz∂uz⎥⎢⎥∂x∂y∂z⎣⎦这里应变张量e是对称的(eij=eji),可表达为:
⎡∂ux⎢
∂x⎢
⎢1∂uy∂uxe=⎢(+)
1∂ux∂uy
(+)2∂y∂x
∂uy
1∂ux∂uz⎤
(+)⎥2∂z∂x⎥1∂uy∂uz⎥(+)⎥ (2.12) ⎢2∂x
∂y∂y⎢⎢1(∂uz+∂ux1∂uz∂uy⎢⎣2∂x
∂z)
2(∂y+∂z
)这里旋转张量Ω是反对称的,可表达为:
⎡⎢
01∂ux∂uy
⎢
2(∂y-∂x)Ω=⎢⎢1∂u⎢
-x∂uy
2(∂y-∂x)0
⎢⎢-1(∂ux-∂uz1∂u
∂u⎢2∂z
)-(y-z)
⎣∂x2∂z∂y读者可对e+Ω=J作检验。
2.2.2物理解释
2∂z∂y⎥
∂uz
⎥∂z⎥
⎥⎦
1∂ux∂uz2(∂z-∂x)⎤⎥⎥12(∂uy∂z-∂uz∂y)⎥⎥ ⎥
⎥0
⎥⎥⎦
2.13) (
为简单起见,考虑平面的几何问题来讨论形变分量和位移之间的关系。如图,经过弹性体内的任意一点P,沿x轴和y轴的正方向取两个微小长度的线段PA=dx,PB=dy。假定弹性体受力后,P,A,B分别移到P’,A’,B’。首先考察线段PA
∂u
和PB的线应变,设P在x方向的位移为u,则A点在x方向的位移将是u+dx,
∂x在y方向的位移为v+
∂v∂udx;B点在x方向的位移将是u+dy,在y方向的位移∂x∂y
为v+
∂u∂u∂v
dy。则PA线段的线应变为其伸长增加量dx与其长度dx之比,为,
∂x∂x∂y
同理,PB线段的线应变为
∂v
。 ∂y
对于上述二维情况,若以PA,PB看做是矩形的两个边,其面积为dxdy,则变形后的P’A’
的长度为dx+
∂u∂v
dx(忽略其他高阶小量),P’B’的长度为dy+dy(忽略其他高阶小量),∂x∂y
则小体积元的面积趋于零时,面积增量与原来面积的极限为:
∂u⎫⎛∂v⎫⎛dx+dxdy+dy⎪-dxdy ⎪
∂x⎭ ∂yS-S0⎝∂u∂v⎝⎭ =≈+
S0dxdy∂x∂y
同样的,对于三维情况,可以得到体积增量和原来体积之比的极限为:
∂u⎫⎛⎛dx+d⎪x
∂x⎭ V-V0⎝⎝=
V0
∂v⎫⎛∂w⎫d+dy+⎪-dz ∂y⎪∂z⎭⎭⎝dxdydz
dxdydz
∂u∂v∂w≈++=divu ∂x∂y∂z
因此(2-12)式对角线元素之和为∆=
∂ux∂uy∂uz
++=tr[e]=∇⋅u为体应变。 ∂x∂y∂z
在线性代数中表示为tr[e]=e11+e22+e33是e的迹。注意,膨胀是由位移场的散度给出的。
现在考察线段PA与PB之间的直角的改变,该改变量有两部分贡献:一部分为由y方向的位移v引起的,即x方向的线段PA的转角α,另一部分为由x方向的位移u引起的,即y方向的线段PB的转角β.由图的关系不难看出:
α≈tgα=
∂v∂u,β≈tgβ= 。在三维情况下也具有同样道理。因此e矩阵中的对角线元∂x∂y
素为线应变,非对角线元素为线段偏转角度之和的一半(采用弧度表示)。
Fig 3
参考上图,一个质点P﹙y、z﹚向逆时针方向扭转到P’﹙y、z﹚,扭转角度为ωx,若其扭转半径为r,根据几何关系可得到:
y
y=rcosα,z=rsinα
其位移形变为
v=-rωxsinα=-zωx w=rωxcosα=yωx
将其分別对y及z微分且相加,得出
ωx=
1⎛∂w∂v⎫
-⎪
2⎝∂y∂z⎭
1⎛∂u∂w⎫1⎛∂v∂u⎫
--⎪,这正是(2-13)式的几个元素 ⎪和ωz= 2⎝∂z∂x⎭2⎝∂x∂y⎭
同理得到ωy=
另外,我们在高等数学中学过旋度的表示为
i
rotu=∇⨯u=
∂∂xux
j∂∂yuy
k∂∂zuz
⎛∂u∂uy⎫⎛∂ux∂uz⎫⎛∂uy∂ux⎫= z---⎪i+ ⎪k⎪j+
∂y∂z∂z∂x∂x∂y⎭⎝⎝⎭⎝⎭
⎛∂u∂uy⎫⎛∂ux∂uz
因此Ω表示了作为刚体的旋转,其中 z-为以x轴为轴的旋转量,-⎪ ∂z⎭∂x⎝∂z⎝∂y
⎫
⎪为⎭
⎛∂uy∂ux⎫
-以y轴为轴的旋转量, ⎪为以z轴为轴的旋转量。因此(2-13)就代表刚度的旋⎝∂x∂y⎭
转。
下面以二维应变为例列举几种应变状态,以此可以理解不同应变状态下的变形。
图
二维情况下不同应变状态的变形举例
2.3 本征值和本征向量
在上式可以看到,像应力张量那样,应变张量是对称的,包含6个独立的参
ˆ,可以找到应变主轴,即:
数。通过计算位移的取向n
ˆ=enˆ
(2.18) u=λn
这里的意义可以仿照应力张量的形式自然得知,在此不再赘述。
在下节我们将发现用脚标符号有助于表达应变张量。方程(2.12)可写成:
1
eij=(
∂iuj+∂jui) (2.20)
2这里假定i和j的顺序为1至3(按x,y,z的取向)。我们用了符号∂xuy=
2.2.1 应变值
∂uy∂x
。
由于应变表示长度的变化除以长度,故应变是无量纲的。与远场地震波通过相联系的动态应变,基本上小于10-6。
2.3 线性的应力—应变关系
在弹性介质中,应力和应变通过应力—应变的本构关系联系起来。应力张量和应变张量之间最一般的线性关系可以写成:
τij=cijklekl=
k=1,3l=1,3
∑∑c
ijkekl
e (2.21)
这里cijkl叫做弹性张量。在此,我们开始采用脚标符号求和的惯用做法。在乘积中任何重复的脚标都意味着求和是从脚标1到3。方程 (2.21)假定是完全弹性的,当应力作用,材料发生变形时,能量没有损失和衰减(有时,这些效应可由。在第六章之前,我们不考虑非弹性状态和衰减。 cijkl为复数来模拟)
应力张量是一个有81个(34)元素的四阶张量。
然而,由于应力张量和应变张量的对称性,它们各自只有六个分量,此时,Cijkl = Cjikl = Cijlk = Cjilk。因此如果不讨论预应力问题,应力分量和应变分量之间的表达式可以表示为:
⎧τ11=c1111e11+c1122e22+c1133e33+c1112e12+c1123e23+c1113e13⎪τ=ce+ce+ce+ce+ce+ce
[***********][***********]⎪22
⎪⎪τ33=c3311e11+c3322e22+c3333e33+c3312e12+c3323e23+c3313e13
(2.22) ⎨
τ=ce+ce+ce+ce+ce+ce[***********][***********]⎪12
⎪τ23=c2311e11+c2322e22+c2333e33+c2312e12+c2323e23+c2313e13⎪⎪⎩τ13=c1311e11+c1322e22+c1333e33+c1312e12+c1323e23+c1313e13
上式只有36个弹性常数。可以证明,对于一个保守系统(即无能量损失),
C1122=C2211,C1133=C3311,C2233=C3322,C1213=C1312,C1223=C2312,C1323=C2313,C1121=C2111,C1113=C1311,C1123=C2311,C2212=C1222,C2213=C1322,C2223=C2322
,
因此对于极端的各向异性介质,C3321=C2133,C3313=C1333,C3323=C2333 (证明从略)。
独立的弹性参数为21个,这些元素只有21个是独立的。这21个元素是确定弹
性固体的最一般形式的应力—应变关系所必须的。
这样一种固体的性质可能随方向变化,如果是这样的话,就称这种介质是各向异性的。与此相反,各向同性的固体,所有方向的性质是相同的。对地球内部,大多数情况下,各向同性是合适的一级近似。但在一些地区观测到各向异性,这是现代研究的一个重要领域。
如果我们作了各向同性的假定,独立参数减至2个:λ和μ,这些参数又叫做介质的拉梅参数。
Cijkl=λδijδkl+μ(δikδjl+δilδjk)
可
见
,
C1
1
=
1
C1=C=μ2λ222
,+
C1122=C1133=C2233=λ
,,,
C1313=C1212=C2323=μ
C1112=C1113=C1123=C2212=C2213=C2223=C3312=C3313=C3323=0C1213=C1223=C1323=0。
我们将看到,拉梅参数及介质的密度将最终确定介质的地震波速度。对各向同性固体,应力—应变方程(2.21)为:
τij=⎡⎣λδijδkl+μ(δilδjk+δikδjl)⎤⎦ekl=λδijekk+2μeij (2.23)
这里我们用eij=eji合并μ的项。注意ekk=tr[e]是e的对角线元素的和。用这方程,我们根据应变可直接写出应力张量:
⎡λtr[e]+2μe112μe12⎢
τ=⎢2μe21λtr[e]+2μe22
⎢2μe2μe32⎣31
⎤⎥
2μe23⎥ (2.24) λtr[e]+2μe33⎥⎦2μe13
μ叫做剪这两个拉梅参数完全描述了各向同性固体里的线性的应力—应变关系。切模量,是介质抗剪切的度量。它的值是所作用的剪切应力与所导致的剪应变比率的一半,即μ=τxy/2exy(这里的2与应变张量的1/2有关)。另一个拉梅参数
λ,没有简单的物理解释。各向同性固体的其他普遍使用的弹性参数中很多具有
特定的物理意义。
2.4 弹性参数的物理意义及相互转换
杨氏模量E:两端拉伸的柱体的张应力与张应变的比率,
正应力和线应变之间满足胡克定律:
σε=
E
所以E表示弹性材料抵抗拉张(或压缩)的能力,E越大说明弹性体越难变形。 Possion比
弹性体纵向伸长(或缩短)∆L后,其横向也会有∆d的变化,
∆d
γ=-
叫做Possion系数。式中负号表示纵向应变和横向应变方向相反,为保证γ为正而取的。Possion系数仅与材料本身的性质有关。实验表明,对于一切介质,γ介于0到0.5之间,对于地幔介质,常用0.25表示地幔的大部分,对于地球外核(液态),取为0.5.
体积模量K:在实际地球中,只受单向压力或拉力的情形很少,一般情况下是各个方向都受力,最常见的是液体静压力。弹性体在静压力P作用下提及变化∆V,则有关系
∆V
P=-K
V
因此体变模量K为流体静压力与由其所导致的体积变化的比率,是介质不可压缩性的度量。
应指出,体变模量K可以从杨氏模量E和Possion系数推导出来。设一个立方体边长为L,体积则为V=L3,在液体静压应力下,整个体积的相对变化率可由
∆V3∆L
=边长的相对变化率决定,即:,就任何一个方向而言,在这个方向上VL
的压应力,使其长度改变∆L'而在与此方向垂直的另两个方向的压应力,又使其长度改变了-2γ∆L',故总的效果是这三个方向压应力作用之和:∆L=(1-2γ)∆L',从而得出:
∆V3∆L∆L'
==3(1-2γ),注意,∆L'为仅在单向压力作用下的长度的改VLL
变量,满足关系式:
∆L'P∆L'=,P=E LEL这样可得到:
ε=
K=
E
31-2γ弹性系数之间的关系: 单向拉伸实验:
τ11=σ,τ22=τ33=τ23=τ13=τ12=0
E
由胡克定律可得: e11=
τ11
,e22=e33=-γe11=-γ
τ11
E
τ11=2μe11+λekk0=2μe22+λekk 0=2μe33+λekk
三式相加得:
ekk=τ11 3λ+2μμ(3λ+2μ) λ+μ代入上式第一式得:E=
代入第二式或第三式得:γ=
各向均匀压缩试验: λ 2(λ+μ)
τ11=τ22=τ33=-p,τ12=τ13=τ23=0 体变模量定义为:K=
由广义胡克定律: -p ekk
-p=2μe11+λekk
-p=2μe22+λekk
-p=2μe33+λekk
三式相加得:
-3p=(3λ+2μ)ekk 与体变模量的定义比较可得K=λ+
在地震学中,我们主要涉及压缩(P)和剪切(S)波的速度。如后面第三章所指出,它们可由这些弹性常数和密度ρ来计算:
P波的速度α表达为: 2μ (2.26) 3
α=
S波的速度β表达为: +2 (2.28) ρ
β= (2.29) ρ
泊松比往往用来度量P波和S波速度的相对大小,表达为:
α2-2β2
(2.30) γ=222(α-β)
注意γ无量纲,在0-0.5之间变化,下限相当于流体的情况(μ=0)。对泊松固体,λ=μ,σ=0.25,α/β=3。多数地壳岩石,泊松比在0.25与0.30之间。
2.3.1 弹性模量的单位
拉梅参数、扬氏模量、体积模量都有像应力那样的单位(帕斯卡)。回顾一下:
1帕斯卡(pa)=1牛顿米-2=1千克米-1秒-2
注意当其被密度除时,结果即为速度平方(适合于方程2.28和2.29)的单位。
练习:
0-1⎫⎛1 ⎪1、地壳中某点的应力状态在北东下坐标系中表示为:τij= 0-10⎪,
-101⎪⎝⎭
单位MPa。求其三个主应力及其方向。
2、已知地壳中某点的主应力为σ1=75bar, σ2=50bar, σ3=-50bar,一斜截面的法线与三个主轴成等角,求该面上的正应力和剪应力。(提示该面的方向余弦为1
,1
,1
3)。
⎛111⎫ ⎪3、已知某点的应力张量为:τ= 111⎪,求其主应力大小和主轴方向。
111⎪⎝⎭
2.1 用方程(2.4),(2.18)和(2.24)说明对各向同性介质,应力主轴总是恰好与应变主轴重合。
2.2 由方程(2.28)和(2.29),根据地震波速度和密度导出拉梅参数的表达式。
2.3对于均匀各向同性弹性体,用应变分量表示应力分量的胡克定律τij=λδijekk+2μeij,试利用弹性常数之间的关系导出应力分量表示为应变分量的表达式:eij=-γ
Eδijτkk+1+γτij。 E
2.3 对泊松比为0.30的岩石,P/S波的速度比是多少?
2.4 观测到实验室里的花岗岩样品的P波速度为5.5km/s,密度为2.6Mg/m3。假定样品是泊松固体,求拉梅参数,扬氏模量、体积模量的值。给出你的以帕斯
卡为单位的答案。如果样品以地球24公里深度所存在的压力压缩,那么这样品体积的相对变化是多少?对这问题,假定当样品压缩时,体积模量没有变化。
2.5 用PREM模型的值(附录1)计算(a)核幔边界(CMB)和(b)内核边界(ICB)两边的体积模量的值。以帕斯卡单位表达你的答案。
2.6 圣地亚哥,加利福尼亚大学在圣地亚哥东北山脉(在Anza附近)建立了Pinon Flat观测台(PFO),仪器包括测量地壳变形的高质量的应变计。
(a) 假定在PFO底下5公里,地震波的速度α=6km/s,β=3.5km/s,密度ρ=2.7Mg/m3,根据这些参数计算拉梅参数λ,μ的值。以帕斯卡单位给出答案。
(b)1922年南加利福尼亚PFO以北80公里的Landers地震(Ms=7.3)后,PFO应变仪观测到相对于地震前,应变有一个大的静态变化。应变张量的水平分量按以下数量变化:e11=-0.26⨯10-6,e22=0.92⨯10-6,e12=-0.69⨯10-6,这里脚标1表示东,2表示北,拉张为正。你可以假定应变的变化是在地震时瞬间出现的。假定在深部,这些值也是正确的,用你在(a)所得到的结果来确定由于地震,在5公里深的地方造成的应力变化,即计算τ11,τ22和τ12的变化。把这看成是假定在垂直方向应变没有变化,应变与深度无关的二维问题。
(c)计算在PFO因Lander地震,应变主轴的取向(水平)。用方位角表达你的答案(北东多少度)。
(d)在PFO观测到应变稳定的长期变化,一年的变化量为:
注意到长期的应变变化接近e11=0.101⨯10-6,e22=-0.02⨯10-6,e12=0.005⨯10-6。
于简单的E-W向拉伸,假定在过去100年,应变的速率是稳定的,初始应力为零,计算5公里深度应力张量的分量(注意,这或许是不很真实的假设!)不包括在5公里深度应力的流体静压分量。
图2.4 在1992年南加利福尼亚地震(MS 7.3)震中
以南80公里的PFO观测台观测到的应变变化。
(e)农民Bob在PFO附近有一平方公里的一小片土地,他将其围起来,并作高精度的测量。农民Bob每年增加或损失多少土地?由于Landers地震,他增加或损失多少土地?
(f)(计算)编写计算机程序计算垂直断层两侧从0°到170°之间不同方位(从北到东以10°增加)的应力。对你在(b)和(d)中计算的应力张量,作表列出断层方位和相应的剪应力及断层两边的法应力(对Landers地震,这些值是应力的变化,不是绝对值)。对每种情况,什么方位有最大的剪应力?
对于(b)中的应力状态:
Vp=6000;Vs=3500;dens=2700;
miu=dens*Vs^2;
lam=dens*Vp^2-2*miu;
e11=0.26e-6;e22=0.92e-6;e12=-0.69e-6;
ekk=e11+e22;
s11=lam*(ekk)+2*miu*e11;
s12=2*miu*e12;
s22=lam*(ekk)+2*miu*e22;
s=[s11,s12;s12,s22];
for theta=0:10:170
N=[cos(deg2rad(theta)) -sin(deg2rad(theta));sin(deg2rad(theta))
cos((deg2rad(theta)))];
theta
S=N*s*N'
end
对于(d)的情况,改变应变值即可。
(g)(计算)一些研究(例如,Stein等人,1992,1994;Harris和Simpson,1992;Harris等人,1995)通过假定沿断层地震破裂的可能性与库伦破裂函数CFF有关,模拟了大地震的空间分布。忽略了孔隙流体压力效应,CFF的变化可以表达为:
∆CFF=∆s+μ∆τn
这里τs是剪应力(曳引力),τn是法应力,μ是静摩擦系数(不能同剪切模量混淆)。注意当剪应力增大,断层上的压应力减小时,CFF增大(回顾一下,在我们的符号约定中,张应力为正,压应力为负)。假定μ=0.2,修改你的计算机程序来计算每个断层取向的∆CFF。作出由于每个断层方位长期的变化,∆CFF年变化的表。
Vp=6000;Vs=3500;dens=2700;
miu=dens*Vs^2;
lam=dens*Vp^2-2*miu;
e11=0.26e-6;e22=0.92e-6;e12=-0.69e-6;
ekk=e11+e22;
s11=lam*(ekk)+2*miu*e11;
s12=2*miu*e12;
s22=lam*(ekk)+2*miu*e22;
s=[s11,s12;s12,s22];
for theta=0:10:170
N=[sin(deg2rad(theta)) cos(deg2rad(theta))];
t=s*N';
sigmaN=N*t;
shear=sqrt(sum(t.^2)-sigmaN^2);
theta
CFF=abs(shear)+0.2*sigmaN
end
(h)(计算)现在假定当他们的长期库伦破裂函数CFF达到某临界阀值时,断层将破裂。至下次地震的时间变量可表达为:
∆t=CFF1000+L-CFF1000 CFFa
这里CFFa是CFF的周年变化,CFF1000是CFF的千年变化,CFF1000+L是CFF的千年变化+Landers地震的变化(注意CFF1000+L≠CFF1000+CFFL)。用在下次地震前每年断层取向加速或减速的时间计算Landers地震的效应。用在下次地震前,加速为正的时间,减速为负的时间的符号约定,陈述你的每年的答案(敬告:这是复杂的)。把断层面上剪应力值的符号反过来,审核你的答案。通常(但不总是)当长期的和Landers地震的剪切变化符号一致(两者都是正或都是负)时,地震时间将提前,当剪切变化的符号不一致时,地震时间将推迟。
(i)提示:对(f)-(h),尤其是(h)作符号校正可能是复杂的。应力可能是正或负。为了得到正确的结果,对每个断层方位规定两个单位矢量,一个平
ˆ)ˆ计算牵引力矢量。行断层(fˆ),另一个垂直于断层(p。通过把应力张量乘以p
ˆ的点积把牵引力矢量分解为剪应力和法应力。然后通过分别计算fˆ与p自然,此
ˆ应是单位长度。 时fˆ和p
(j)Landers地震后,没有观测到PFO附近地震活动性(小震活动)增强。这说明临界值的CFF模型是正确的吗?