2 第二章应力和应变_wan

1、 第二章应力和应变地震波传播的任何定量的描述，都要求其能表述固体介质的内力和变形的特征。现在我们对后面几章所需要的应力、应变理论的有关部分作简要的复习。虽然我们把这章作为独立的分析，但不对许多方程进行推导，读者想进一步了解其细节，可查阅连续介质力学的教科书。三维介质的变形称为应变，介质不同部分之间的内力称为应力。应力和应变不是独立存在的，它们通过描述弹性固体性质的本构关系相联系。应力的表述应力张量2.1.1应力表示考虑一个在静力平衡状态下，均匀弹性介质里一个任意取向的无限小平面。平面的取向可以用这个平面的单位法向矢量n来规定。在n方向的一侧施加在此面单位面积上的力叫做牵引力，用矢量t二(t,t

2、,t)表示。在n相反方向的另一侧施加在此面上的xyz力与其大小相等，方向相反，即t(-n)=-t(n)。t在垂直于平面方向的分量叫做法应力，平行于平面方向的分量叫做剪应力。在流体的情况下，没有剪应力，t=-pn，这里P是压强。在笛卡尔坐标系（图2.1）里，应力张量T可以用作用于yz,xz,xy平面的牵引力来定义+：t(x)t(y)tTTTXxxxxxyxzt(x)t(y)t=TTTyyzyxyyyzt(x)t(y)tTTTzzz1zxzyzzT二2.1)应力分量的符号规定如下：对于正应力，我们规定拉应力为正，压应力为负。对于剪应力，如果截面的外法线方向与坐标轴一致，则沿着坐标轴的正方向为正，反

3、之为负；如果截面方向与外法线方向相反，则沿着坐标轴反方向为正。图2.1在笛卡尔坐标系里描述作用在无限小立方体面上的力的牵引力矢量t（X）,t（y）,t（z）。+通常，应力张量用脚标互换的（2.1）式来定义，第一个脚标表示面的法向方向。实际上，当工是对称的时，没有什么差别。因为固体处于静力平衡状态，所以没有因剪应力作用而产生的净扭转。例如，当考虑XZ平面的剪应力时，为了使扭转相互抵消，必须有P=T。类似地有:xzzxP=P,P=P。故应力张量是对称的，即：2.2)xyyxyzzyPPPxxxyxzPPPxyyyyzPPPxzyzzz应力张量只包含6个独立的要素，它们足以完全描述介质中一个给定点的

4、应力状态。2、任意一个面上的应力可以由应力张量表示一点的九个应力分量如果能够完全确定一点的应力状态，则其必须能够表达通过该点的任意斜截面上的应力矢量。为了说明这一问题，在O点用三个坐标面和一任意斜截面截取一个微分四面体单元，斜截面的法线方向矢量为n，它的三个方向余弦分别为l，m和n。设斜截面上的应力为pn，i，j和k分别为三个坐标轴方向的单位矢量，卩在坐标轴上的投影分别为p,n.xp,p。则应力矢量可以表示为yzpn=pxi+pyj+pzk同样，把单位体积的质量所受的体积力Fb沿坐标轴分解，有Fb=Fbxi+Fbyj+Fbzk设S为AABC的面积，贝UNOBC=lS,AOCA=mS,AOAB=

5、nSAABC的法线方向的单位矢量可表示为n=li+lj+mk微分四面体在应力矢量和体积力作用下应满足平衡条件，设h为O点至斜面ABC的高，由x方向的平衡，可得近=0p-BC-zAOAC-tAOAB+=0将公式AOEC二屁4OCA二漱,心E二恋代入上式，则p.+二对于微分四面体单元，h与单元体棱边相关，因此与1相比为小量，趋近于零（在忽略体力的情况下，F=0），因此bx=4+叨+7同理Py二坯+6严+厂/Ps二%丿+匸即腕+耳充如果采用张量记号，则上述公式可以表示为代=疔化上式给出了物体内一点的9个应力分量和通过同一点的各个微分面上的应力之间的关系。这一关系式表明，只要有了应力分量，就能够确定一

6、点任意截面的应力矢量，或者正应力和切应力。因此应力分量可以确定一点的应力状态。任一个取向由n定义的平面，一个侧面作用于另一个侧面的牵引力为应力张量与n的乘积，即：t(n)=Tn=t(n)TTTnxxxxyxzxt(n)=TTTnyxyyyyzyt(n)TTTnzxzyzzzz2.3)这可以通过对由垂直于n的平面和xy,xz,yz平面所围限的四面体（柯西四面体）面上的力求和作出说明。简单地说，应力张量是给出相对于法向矢量n的牵引力矢量t的线性算子,从这个意义上来说，应力张量与任何特定的坐标系无关。在地震学中，我们几乎总是把应力张量写成笛卡尔几何学里的一个3x3的矩阵。注意到对称的要求，应力张量的

7、独立参数由9个减少为6个，呈现为最一般形式的二阶张量（标量为零阶张量，矢量为一阶张量，等等）。应力张量通常随在物质里的位置而变化，它是作用在固体里每一点的无限小的面上的力的度量。应力只给出了这些面由一边作用于另一边的力的度量，计量标准是单位面积上的力。然而，可能有其他力（如重力）。这些力称为体力，计量标准是每单位体积或单位质量上的力。2.1.3坐标变换一点的应力分量不仅随着变形体中点的位置在改变，而且即使在同一点，由于截面的法线方向不同，截面上的应力也不一样。假设已知在坐标系Oxyz中，弹性体中的某点的应力状态表示为：TTTxxxyxzT=TTTxyyyyzTTTxzyzzz则对于新的坐标系O

8、xxyzxlmn1i1ylmn222zlmn333yz，这点的各应力分量是多少就是本小节讨论的问题。设新旧坐标系的方向余弦表如下：根据（2.3）可知，沿十的应力可以表示为:t(x)=Tx=TTTxxxyxzTTTxyyyyzTTTxzyzzzl1m1n1此应力有三个分量，将其再投影到r方向即得到在新坐标系中的各分量:tttt=it(xr)=bmnxxtxytxzt1mxx111xyyyyz1tttnxzyzzz1-c(n1=12T+m2T+n2T+2lnt+2lnt+2mnt1xx1yy1zz11xy11xz11yzttt7一t=jt)=tmnxxtxytxztimxy222xyyyyzitt

9、tnxzyzzz1=11t+mmt+nnt+21nt+21nt+2mnt12xx12yy12zz11xy11xz11yzTmnttt_111一11ixxxyxz123t=1mntttmmmij222xyyyyz1231mntttnnn333xzyzzz123所以=NTtN由此可知一个斜面上的正应力可以表示为：gG)=kn1Ftttnn1xxItxytxztxnxyzxyyyyzytttnxzyzzzz其剪应力为：2.1.4主应力和应力主轴对任何应力张量，总是可以找到一个方向n，使得在垂直于n的面上，没有剪应力，也就是说，t（n）沿n方向，在这种情况下：t（n）=Xri=tntn-九n=o（2.

10、4）（t-i九）n=o这里i是单位矩阵，久是标量。这是一个本征值问题，只有detT-U=0（2.5）才有非无效的解。这是一个有三个本征值九、九、九解的三方程（不要把这些值123同后面将讨论的拉梅参数相混淆）。由于T是对称的，是实数，所以本征值也是实数。相应于本征值的本征矢量n,n,n是正交的。它们定义了应力主轴。垂直于应力主轴的平面叫做主平面。我们通过相似性变换，把t旋转到n,n,席的坐标系里：t1tr=NttN=000t020t32.6)这里TR是旋转的应力张量，T、T、12t是主应力(与本征值九，九，九相等)，N是3123本征矢量矩阵：n(1)n(2)n(3)xxxn(1)n(2)n(3)

11、yyyn(1)n(2)n(3)zzzNt=N-1为归一化到单位长度的正交的本征矢量。2.7)在MATLAB中矩阵的本征值和本征向量的求法为：X,D=eig(sigma);其中sigma为三维应力矩阵，X为三个本征向量，D为本征值矩阵。得到本征向量和本征值后，你可以采用X*sigma*X得到的结果验证是否是向量矩阵。如果T二T二T，那么应力场处于流体静压状态，没有任何取向的面有剪应123力。在流体的情况下，应力张量可写成：-P00T=0-P0这里P是压强。对于垂直应力不变的应力状态,其水平方向不同的应力状态可以表示为TxxTxy0-Txx-cTxy0_TT0，其主应力方程可以表示为TT-c0 x

12、yyyxyyy00T00T-c=0，或者展开为T=zz解这个方程可得到三个主应力为T-爲2-(+T(T-T210,2.8)zzxxyy1Q=121Q=22xxyyxyT+Txxyy+T+Txxyy0及G二T3zz2.1.1应力值应力以单位面积上的力为计量单位，在国际单位制（SI）中的单位是：1帕斯卡（pa）=1牛顿米-2回顾一下，1牛顿=1千克米秒-2=105达因。另一个普遍采用的力的单位是巴：1巴（bar）二105pa1千巴二108pa二100Mpa1兆巴（Mbar）二1011pa二100Gpa（千兆帕）如表2.1用参考模型PREM（Dziewonski和Anderson，1981）所给出的

13、值所示，在地球里压力随深度快速增大。在400公里的深度，压力达到13.4Gpa，在核幔边界达136Gpa，在内核边界达329Gpa。与此对比，月球中心的压力仅为4.8Gpa,相当于在地球150公里深度所达到的值（Latham等人，1969）,这是由于月球的质量小得多。表2.1地球内部压力与深度的关系深度（公里）区域压力（Gpa）0-24地壳0-0.624-400上地幔0.6-13.4400-670过渡区13.4-23.8670-2891下地幔23.8-135.82891-5150夕卜核135.8-328.95150-6371内核328.9-363.9这些压力是地球内部的流体静压力。深部的剪切应

14、力值小得多，包括与地幔对流有关的应力和由地震波传播所产生的动态应力。静态应力可能存在于地壳上部脆性部分。测量地壳剪应力是现代研究的课题。应力的量级是有争议的问题，地壳应力可能在100巴和1000巴（10-100Mpa）之间，在接近活动断层的区域,应力有降低的趋势（活动断层的作用降低了应力）。2.2应变张量现在让我们来考虑怎样描述连续介质里点的位置的变化。任何一点与参考时间t的相对位置都可以用一个矢量场来描述，即位移场为：0u（r）二rr（2.9）00这里r是现在的位置，r是参考点的位置。位移场是一个重要的概念，在这本书0中经常涉及到。它是位置变化的绝对度量。与此相对照，应变是位移场相对变化的局

15、部度量，即位移场空间梯度的度量。应变与材料的变形或形状的变化有关，而与位置的绝对变化无关。例如，张应变是按照长度的相对变化来定义的。如果把一根100米长的细绳，一端固定，在另一端均匀地拉长到101米，那么沿细绳位移场从0变化到1米。依此，在绳的任何地方，应变场为0.01(1%)的常数。考虑离开参考位置x一个小的距离的某一点x的位移u二(u,u,u)，对u0 xyz作泰勒级数展开，得到：dududuu(x)=dxx-dyxdzdududu贰dy牙dududu=u(x)+0dxd=u(x)+Jdy0dz2.10)这里d=x-x0，假定偏导数土,土等很小，它们的乘积可以忽略。那么，根据无限小应变理论

16、，我们就可以忽略展开式中的高阶项。庆幸的是在地震学中地球的应变几乎总是很小的，以致于这种近似是恰当的。我们可以通过把J分成对称和反对称部分，把刚性旋转部分离出来：dududux-dxx-dyx-dzdududuJ=A=e+Qdxdydzdudududxdyzdz2.11)这里应变张量e是对称的(e=e)，可表达为：ijjiduxdxdudu、e=(+)dxdy HYPERLINK l bookmark62 dudu(壬+y) HYPERLINK l bookmark64 dydxdu1dudu卞)1du2(苛+du)dy2.12)1dudu2(盂+dZ)du(zdydu+牙)duzdz这里旋转

17、张量Q是反对称的，可表达为:(dudyxduydxdudu(xz)dzdxduQu(x-斗)dydx1du2(1Zduy2.13) HYPERLINK l bookmark78 dudu(x于)dzdx1du2(1Zduy读者可对e+Q=J作检验。可通过考虑一个无限小的立方体所发生的情况（图2.2）对e和Q效应作说明。e的非对角线元素形成了剪应变。例如，在二维的情况下，如果0二0和没有体积变化，那么有dududuxz0,xdxdzdzdu,以及Je0duzdxduxdz02.14)图2.2用x-z平面的正方形的图解说明应变张量e和旋转张量0的不同效应。e的非对角线元素形成了剪应变（左边的正方形

18、），而0形成了刚性旋转（右边的正方的）。相对于有效无限小应变理论，这里所标的变形是大大放大了的。这里9是每边旋转的角度（弧度，不是度数！），注意边之间角度的总变化是29。与此相对照，矩阵0形成刚性旋转。例如，如果e-0，那么有：生七，以及:0-90dudXduxdz02.15)在这两种情况下，材料没有体积变化。相对的体积增大，或膨A=（V-V）/V00由x，y，z方向的拉张之和给出：TOC o 1-5 h z HYPERLINK l bookmark70 dududu HYPERLINK l bookmark130 A=x+a+ztre=V-udxdydz这里tree+e+e是e的迹。注意，膨

19、胀是由位移场的散度给出的。112233那么，位移场的旋度是什么呢？我们回顾一下矢量场的旋度的定义 HYPERLINK l bookmark132 dudududududu HYPERLINK l bookmark134 Vxu(土y)X+(xw)y+(ax HYPERLINK l bookmark136 dydzdzdxdxdy2.16)2.17)把方程与（2.13）比较说明，只要Q不为零，那么Vxu就不为零，位移场就包含某种刚性旋转。像应力张量那样，应变张量是对称的，包含6个独立的参数。通过计算位移的取向n，可以找到应变主轴，即：u九neni（2.18）图2.3样品沿x方向的张应变导致剪应变

20、。内部的角度是变化的。这与上节讨论应力张量的情况类似。三个本征值为主应变e，e，e，而本征矢量123规定了应变主轴。注意，除了eee的流体静应变的情况外，总是有某种剪123应变出现。例如，考虑一个仅在x方向拉张的正方形（图2.3），那么e可表达为：Que0 x-01=Qx00|_00e=2.19)平行于坐标轴的线之间的角度没有变化，但中间角度的线发生旋转。如果把坐标系旋转45，那么与剪应变相联系的角度的变化是很明显的（在这种情况下，e有两个非对角线项）。在下节我们将发现用脚标符号有助于表达应变张量。方程（2.12）可写成：1e=（du+du）（2.20）ij2ijjiQu这里假定i和j的顺序为

21、1至3（按x,y,z的取向）。我们用了符号Qu=亍。xyQx2.2.1应变值由于应变表示长度的变化除以长度，故应变是无量纲的。与远场地震波通过相联系的动态应变，基本上小于10-6。2.3线性的应力应变关系在弹性介质中，应力和应变通过应力应变的本构关系联系起来。应力张量和应变张量之间最一般的线性关系可以写成：t二ce二Hce（2.21）ijijklklijkekl这里c叫做弹性张量。在此，我们开始采用3脚标符号求和的惯用做法。在乘积ijkl中任何重复的脚标都意味着求和是从脚标1到3。方程（2.21）假定是完全弹性的，当应力作用，材料发生变形时，能量没有损失和衰减（有时，这些效应可由c为复数来模拟

22、）。在第六章之前，我们不考虑非弹性状态和衰减。ijkl应力张量是一个有81个（34）元素的四阶张量。然而，由于应力张量和应变张量的对称性，它们各自只有六个分量。因此如果不讨论预应力问题，应力分量和应变分量之间的表达式可以表示为：TOC o 1-5 h zt=ce+ce+ce+ce+ce+ce11111112221333141215231631t=ce+ce+ce+ce+ce+ce211122222333241225232631t=ce+ce+ce+ce+ce+ce33311132223333341235233631（2.22）t=ce+ce+ce+ce+ce+ce124111422243334

23、41245234631t=ce+ce+ce+ce+ce+ce511152225333541255235631t=ce+ce+ce+ce+ce+ceV31611162226333641265236631上式只有36个弹性常数cCj=1,2,3,4,5,6）。可以证明，对于一个保守系统(及无能量损失)c二c(证明从略)。因此对于极端的各向异性介质，独立的弹(/jTOC o 1-5 h z性参数为21个(c=c,c=c,c=c,c=c,c=c、12211331144115511661c二c,c二c,c二c,c二c、c二c,c二c,c二c、c二c,c二c、23322442255226623443355

24、3366345544664c二c)。，这些元素只有21个是独立的。这21个元素是确定弹性固体的最一6556般形式的应力应变关系所必须的。这样一种固体的性质可能随方向变化，如果是这样的话，就称这种介质是各向异性的。与此相反，各向同性的固体，所有方向的性质是相同的。对地球内部，大多数情况下，各向同性是合适的一级近似。但在一些地区观测到各向异性，这是现代研究的一个重要领域(参见11.3节，关于各向异性)。如果我们作了各向同性的假定，此时，独立参数减至2个:九和卩，这些参数又叫做介质的拉梅参数。我们将看到，拉梅参数及介质的密度将最终确定介质的地震波速度。对各向同性固体，应力应变方程(2.21)为：t5

25、+卩(55+55)e二九e+2ye(2.23)i/i/klil/kik/lkli/kki/这里我们用e二e合并卩的项。注意e二tre是e的对角线元素的和。用这方i/ikk程，我们根据应变可直接写出应力张量：九tre+2ye112ye212ye312ye12九tre+2ye222ye322ye132ye23九tre+2ye332.24)这两个拉梅参数完全描述了各向同性固体里的线性的应力应变关系。y叫做剪切模量，是介质抗剪切的度量。它的值是所作用的剪切应力与所导致的剪应变比率的一半，即y=T/2e。另一个拉梅参数九，没有简单的物理解释。各向同xyxy性固体的其他普遍使用的弹性参数中很多具有特定的物

26、理意义。3弹性参数的物理意义及相互转换杨氏模量E：两端拉伸的柱体的张应力与张应变的比率，正应力和线应变之间满足胡克定律：=E所以E表示弹性材料抵抗拉张(或压缩)的能力，E越大说明弹性体越难变形。Possion比弹性体纵向伸长（或缩短）AL后，其横向也会有Ad的变化，y_-d叫做Possion系数。式中负号表示纵ALl向应变和横向应变方向相反，为保证Y为正而取的。Possion系数仅与材料本身的性质有关。实验表明，对于一切介质，Y介于0到0.5之间，对于地幔介质，常用0.25表示地幔的大部分，对于地球外核（液态），取为0.5.体积模量K:在实际地球中，只受单向压力或拉力的情形很少，一般情况下是各

27、个方向都受力，最常见的是液体静压力。弹性体在静压力P作用下提及变化AV，则有关系AVP_-KV因此体变模量K为流体静压力与由其所导致的体积变化的比率，是介质不可压缩性的度量。应指出，体变模量K可以从杨氏模量E和Possion系数推导出来。设一个立方体边长为L,体积则为V=L3,在液体静压应力下，整个体积的相对变化率可由边长的相对变化率决定，即：AV_3AL，就任何一个方向而言，在这个方向上VL的压应力，使其长度改变AL而在与此方向垂直的另两个方向的压应力，又使其长度改变了-2yALf，故总的效果是这三个方向压应力作用之和：AL_6-，从而得出：_3（1-2y）AL，注意，AL为仅在单向压力作用

28、下的长度VLL的改变量，满足关系式：AL_,P_ELEL这样可得到：弹性系数之间的关系单向拉伸实验:T_G,T_T_T_T_T_0112233231312TTe_ii,e_e_ye_71111E223311E由胡克定律可得：T二2pe+九e1111kk0二2pe+九e22kk0二2pe+九e33kk三式相加得：Tekk=&+k）代入上式第一式得：E=呻+2九+卩代入第二式或第三式得：Y=2（2+卩）各向均匀压缩试验：T=T=T=p,T=T=T=0112233121323体变模量定义为：K=Pekk由广义胡克定律：p=2pe+2e11kkp=2pe+2e22kkp=2pe+2e33kk三式相加得

29、：3p=&+2応kk22.26)与体变模量的定义比较可得K=2+2卩在地震学中，我们主要涉及压缩（P）和剪切（S）波的速度。如后面所指出，它们可由这些弹性常数和密度p来计算：P波的速度表达为：2.28)a=；2斗PPS波的速度0表达为:2.29)泊松比往往用来度量P波和S波速度的相对大小，表达为：2.30)a22022(a2-p2)注意丫无量纲，在0-0.5之间变化，下限相当于流体的情况（卩二0）。对泊松固体，九二二0.25,a/Pf：3。多数地壳岩石，泊松比在0.25与0.30之间。2.3.1弹性模量的单位拉梅参数、扬氏模量、体积模量都有像应力那样的单位（帕斯卡）。回顾一下：1帕斯卡（pa）

30、=1牛顿米-2二讦克米-1秒-2注意当其被密度除时，结果即为速度平方（适合于方程2.28和2.29）的单位。练习：10-1、1、地壳中某点的应力状态在北东下坐标系中表示为：工=0-10,ijI-101丿单位MPa。求其三个主应力及其方向。2、已知地壳中某点的主应力为b=75bar,q=50bar,q=-50bar,斜截123面的法线与三个主轴成等角，求该面上的正应力和剪应力。（提示该面的方向余弦为咅吉吉）。(113、已知某点的应力张量为：工=11b11,求其主应力大小和主轴方向。用方程（2.4），（2.18）和（2.24）说明对各向同性介质，应力主轴总是恰好与应变主轴重合。由方程（2.28）和

31、（2.29），根据地震波速度和密度导出拉梅参数的表达式。对于均匀各向同性弹性体，用应变分量表示应力分量的胡克定律T二九6e+2比，试利用弹性常数之间的关系导出应力分量表示为应变分量的ijijkkij表达式：e=-6t+toijEijkkEij2.3对泊松比为0.30的岩石，P/S波的速度比是多少？2.4观测到实验室里的花岗岩样品的P波速度为5.5km/s,密度为2.6Mg/m3。假定样品是泊松固体，求拉梅参数，扬氏模量、体积模量的值。给出你的以帕斯卡为单位的答案。如果样品以地球24公里深度所存在的压力压缩，那么这样品体积的相对变化是多少？对这问题，假定当样品压缩时，体积模量没有变化。2.5用P

32、REM模型的值（附录1）计算（a）核幔边界（CMB）和（b）内核边界（ICB）两边的体积模量的值。以帕斯卡单位表达你的答案。2.6圣地亚哥，加利福尼亚大学在圣地亚哥东北山脉（在Anza附近）建立了PinonFlat观测台（PFO），仪器包括测量地壳变形的高质量的应变计。（a）假定在PFO底下5公里，地震波的速度a=6km/s,卩=3.5km/s，密度p=2.7Mg/m3，根据这些参数计算拉梅参数九,卩的值。以帕斯卡单位给出答案。（b）1922年南加利福尼亚PFO以北80公里的Landers地震（Ms=7.3）后，PFO应变仪观测到相对于地震前，应变有一个大的静态变化。应变张量的水平分量按以下数

33、量变化：e=0.26x10-6,e二0.92x10-6,e=0.69x10-6，112212这里脚标1表示东，2表示北，拉张为正。你可以假定应变的变化是在地震时瞬间出现的。假定在深部，这些值也是正确的，用你在（a）所得到的结果来确定由于地震，在5公里深的地方造成的应力变化，即计算T,1和T的变化。把这112212看成是假定在垂直方向应变没有变化，应变与深度无关的二维问题。（c）计算在PFO因Lander地震，应变主轴的取向（水平）。用方位角表达你的答案（北东多少度）。（d）在PFO观测到应变稳定的长期变化，一年的变化量为：e=0.101x10-6,e=-0.02x10-6,e=0.005x10

34、-6。注意到长期的应变变化接近112212于简单的E-W向拉伸，假定在过去100年，应变的速率是稳定的，初始应力为零，计算5公里深度应力张量的分量（注意，这或许是不很真实的假设！）不包括在5公里深度应力的流体静压分量。3533条地亚哥太平1151611819Landers地震-117洛杉矶图2.4在1992年南加利福尼亚地震(M7.3)震中S以南80公里的PFO观测台观测到的应变变化。农民Bob在PFO附近有一平方公里的一小片土地，他将其围起来，并作高精度的测量。农民Bob每年增加或损失多少土地？由于Landers地震，他增加或损失多少土地？(计算)编写计算机程序计算垂直断层两侧从0到170之

35、间不同方位(从北到东以10增加)的应力。对你在5)和(d)中计算的应力张量，作表列出断层方位和相应的剪应力及断层两边的法应力(对Landers地震，这些值是应力的变化，不是绝对值)。对每种情况，什么方位有最大的剪应力？对于(b)中的应力状态：Vp=6000;Vs=3500;dens=2700;miu=dens*Vs2;lam二dens*Vp2-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;fortheta=0:10:170N=cos(deg2rad(theta)-sin(deg2rad(theta);sin(deg2rad(theta)cos(deg2rad(theta);thetaS=N*s*Nend对于(d)的情况