呂愛鐘， 張曉莉， 王少杰

(華北電力大學 可再生能源學院，北京 102206)

隨著現(xiàn)代科學技術的發(fā)展，復合材料以其重量輕、比強度高、比剛度高、可設計性強等特點，在航空航天、汽車、石油化工、建筑橋梁、艦船以及生物工程等領域得到了越來越廣泛的應用.常用的復合材料板多為各向異性板，由于結構裝配或功能方面的需求，常常需要在板中打孔，而各向異性板在不同方向呈現(xiàn)出不同的力學性能，對于這種開孔結構來說，特別突出的一個問題就是開孔引起的應力集中.在荷載作用下，開孔可能會對結構的承載能力、使用壽命產生嚴重影響，因此，精確計算出孔洞周邊的應力分布對判斷結構的穩(wěn)定性以及采取何種加固措施等具有重要的指導意義.為此，國內外許多學者都對各向異性含孔結構進行了應力分析.

Lekhnitskii[1-2]最早采用復變函數(shù)方法，研究了帶有一個橢圓孔的各向異性板在平面荷載作用下的應力計算問題.在此后很長一段時間內，只有圓孔或橢圓形孔才有精確解.近幾年來，一些學者通過復變函數(shù)方法，也研究了三角形孔[3-6]、矩形孔[6-9]、某六邊形孔[9]和多邊形孔[10]在無窮遠處作用不同荷載時的孔邊應力解析解.Ukadgaonker等[11]對不規(guī)則孔的孔邊應力進行了分析.但是，上述研究成果獲得的解析解只有在孔邊才是適用的[12]，所以上述研究中的算例也只對孔邊應力進行了分析，而沒有涉及孔外域，因此存在一定的局限性.筆者[12]曾對這一領域進行了研究，找到了含任意孔形的正交各向異性隧洞應力計算的精確解析方法，并對馬蹄形和直墻半圓拱形隧洞進行了孔邊應力的分析,但是，并未對各向異性板中含有一些不規(guī)則形狀孔等更復雜的孔形進行應力解析分析，對孔邊應力集中明顯的位置及其孔外域的應力分析也不夠全面.

為了研究各向異性板開孔結構應力分布的特殊性，本文選取具有明顯尖點的某不規(guī)則形狀孔和正六角形孔為例，根據帶孔正交各向異性板的應力解析解，從板不同的纖維方向、開孔形狀和外荷載方向幾個方面分析孔邊及部分孔外域的應力，并且與各向同性開孔結構進行對比，總結正交各向異性板孔邊及孔外域的應力分布規(guī)律.不規(guī)則孔形在文獻[11]中有討論，但因為其解并未涉及孔外域的應力分布，本文探討此孔形一方面是進行對比驗證，更重要的是對孔外域的應力也進行研究，著重分析應力分布的規(guī)律.正六角形孔在工程中應用得也很多，但至今尚未見到其精確的分析結果，其主要原因可能是缺乏可用的高精度的正六角形映射函數(shù).

1 基本方程和應力解析解的求解

工程中常用的復合材料結構在很多場合都具有對稱性，因此在很多情況下，各向異性的問題都可以化為正交各向異性的問題來討論.復變函數(shù)中的保角變換方法是求解孔洞問題的有效工具，對于正交各向異性問題，求解時涉及3套笛卡爾坐標系：z=x+iy，z1=x+μ1y，z2=x+μ2y，通過建立3套極坐標系ζ=ρeiθ，ζ1=ρ1eiθ1，ζ2=ρ2eiθ2可以將z，z1，z23個物理平面的孔外域分別映射到ζ，ζ1，ζ23個像平面的單位圓外域，映射函數(shù)分別用z=ω(ζ)，z1=ω1(ζ1)，z2=ω2(ζ2)表示，其中，z=ω(ζ)可以表示為以下Laurent級數(shù)的形式[13]：

(1)

式中：R和Ck是代表孔的大小和孔形的常數(shù)；k為整數(shù).當k取較大的值時，式(1)可以代表足夠多的孔形.取k的最大值為n.由于3套笛卡爾坐標系z，z1，z2必須滿足如下數(shù)學關系[1-2]：

(2)

式中：γ1=(1-iμ1)/2，γ2=(1-iμ2)/2，δ1=(1+iμ1)/2，δ2=(1+iμ2)/2，μ1=α1+iβ1，μ2=α2+iβ2，αk和βk(k=1,2)是與材料性質有關的材料常數(shù)，且β1>0，β2>0.將式(1)代入式(2)，并假設3套極坐標系在孔邊界所對應的單位圓上滿足ζ=ζ1=ζ2=σ=eiθ，可推導出z1=ω1(ζ1)和z2=ω2(ζ2)的表達式和ζ，ζ1，ζ2之間的關系如下[12]：

(3)

(4)

式(3)對任意形狀的孔形都適用，式(4)對域內(|ζ|≥1)任意一點都成立.

如圖1所示，考慮一個含孔正交各向異性板，承受分布于板邊緣并作用在面內的力，孔形關于整體坐標系的x軸對稱，局部坐標系x′Oy′沿彈性主向.如果板的尺寸遠大于開孔的尺寸，并且開孔不在板邊的附近，那么可以假定板無限大.除此之外，假定作用在板內的荷載遠大于體力，并且無z方向的分量，此時，體力可忽略不計，問題簡化成了無限域中的平面應力問題.根據平衡條件和變形協(xié)調條件可以得到應力函數(shù)F(x，y)所滿足的方程[14]：

(5)

式(5)的解與特征方程(6)的根有關.

a11μ4-2a16μ3+(2a12+a66)μ2-2a26μ+a22=0

(6)

只討論μ1≠μ2的情形，引入2個解析函數(shù)F1(z1)、F2(z2)，則式(5)的解可以表示為

(7)

令Φ1(z1)=dF1(z1)/dz1，Φ2(z2)=dF2(z2)/用Φ1(z1)、Φ2(z2)表示的應力邊界條件為

圖1受面內均布荷載作用且?guī)в腥我庑螤羁椎恼桓飨虍愋园?/p>

Fig.1Anorthotropicplatewithanarbitraryshapedholeunderin-planeloadings

dz2，這樣，求解各向異性板平面應力問題就轉化為尋找滿足相應的邊界條件的應力解析函數(shù)Φ1(z1)和Φ2(z2)的問題.

(8)

式中出現(xiàn)的z1、z2是孔邊點.因為討論的是帶孔無限域問題，且孔邊無外荷載，則f1=f2=0，且Φ1(z1)、Φ2(z2)具有以下形式：

(9)

式中：實常數(shù)B*、B′*、C′*可以通過無限遠處的外荷載(σx∞、σy∞、τxy∞)和與材料性質有關的常數(shù)(α1、α2、β1、β2)求出[15].

(10)

因在孔的邊界有ζ1=σ=eiθ、ζ2=σ=eiθ，所以應力邊界條件式(8)最終可以整理為如下形式：

(11)

根據式(3)、式(10)和式(11)，采用冪級數(shù)解法.將σ-k=coskθ-isinkθ代入式(11)的左右兩邊并展開，可根據sinkθ、coskθ(k=1,∞)前面的系數(shù)相等列出求解ak、bk的線性方程組，具體計算過程略.

所以，一旦正交各向異性板的彈性常數(shù)給定，式(6)中的系數(shù)就可以確定[15]，特征方程的根可解，再根據給定的外荷載，實常數(shù)B*、B′*、C′*也可以求得.結合求解的ak、bk，再根據式(9)和式(10)，解析函數(shù)Φ1(z1)、Φ2(z2)就可以求解.而笛卡爾坐標系下和正交曲線坐標系下的應力分量也可以根據下式進行求解：

(12)

(13)

(14)

式中：φ是局部坐標系x′Oy′與整體坐標系xOy的夾角，見圖1.

至此，對于孔的邊界，由于ρ=1、ζ1=ζ2=ζ=σ=eiθ，應力分量非常容易求出，而對于ζ≠σ的孔外域點，應力求解稍微復雜一些，必須首先根據式(4)由ζ求出對應的ζ1、ζ2，再根據Φ1(z1)、Φ2(z2)進行應力求解.需要注意的是，由式(4)求出的ζ1、ζ2各有2n個解，但由于映射函數(shù)是將物理平面上的孔外域映射到像平面上的單位圓外域，在這2n個解中，只有n個解在單位圓外域.除此之外，再根據物理平面上的2個點映射到像平面上也相鄰的原則，可以找出n個解中符合題意的解.

本文得出的結果是適用于任意形狀孔形和任意纖維角度的精確解，只要孔的解析函數(shù)已知，就可以求出面內荷載作用下孔邊及其周圍的應力場.

2 算例分析和結果比較

2.1 不規(guī)則孔形的孔邊切向應力分析

對于帶有不規(guī)則孔形的正交各向異性板，文獻[11]中考慮了由石墨-環(huán)氧樹脂材料(G-E)制成的具有不同纖維方向([0°/-90.0°]s，[45.0°/-45.0°]s)的復合層壓板，給出的不規(guī)則孔形的映射函數(shù)系數(shù)Ck為：C3=0.06，C6=0.04，C8=0.03，C9=0.03，其他Ck均為零，特征方程的復根為：當[0°/90.0°]s時,有μ1=0+3.640 4i，μ2=0+0.274 7i；當[45.0°/-45.0°]s時,有μ1=-0.859 7+0.510 9i，μ2=0.859 7+0.510 9i.

為了方便進行對比驗證和分析，取上述參數(shù)與文獻[11]完全相同，由于討論的是無限大平板問題，孔的大小對應力場的求解沒有影響，故本文中映射函數(shù)的R值直接取為1.0.當板的無窮遠處分別單獨作用沿x軸或y軸方向的單向均布拉伸荷載時，即圖1中τxy∞=0，而σx∞=σ或σy∞=σ，根據上文介紹的應力解析解的求解方法，通過自編FORTRAN程序，實現(xiàn)該不規(guī)則孔形孔邊的應力分量的計算，結果如圖2和圖3.由于孔邊沒有外荷載作用，所以在孔邊，徑向應力σρ=0，剪應力τρθ=0，故本文圖中虛線繪出的都是孔邊的切向應力與外荷載的比值σθ/σ的分布，其中，θ是ζ平面的極角.符號規(guī)定拉應力為正、壓應力為負.

a 各向同性材料

b G-E[0°/-90.0°]s

c G-E[45.0°/-45.0°]s圖2 帶有不規(guī)則孔形的3種材料結構在σx∞=σ單獨作用下的孔邊切向應力分布

Fig.2Distributionofnormalizedtangentialstressesforthreekindsofplateswithirregularshapedholeunderuniaxialloadingσx∞=σ

a 各向同性材料

b G-E[0°/-90.0°]s

c G-E[45.0°/-45.0°]s圖3 帶有不規(guī)則孔形的3種材料結構在σy∞=σ單獨作用下的孔邊切向應力分布

Fig.3Distributionofnormalizedtangentialstressesforthreekindsofplateswithirregularshapedholeunderuniaxialloadingσy∞=σ

首先將圖2b、2c和圖3b、3c的結果與文獻[11]進行對比驗證，發(fā)現(xiàn)文獻[11]中獲得的孔邊切向應力的分布規(guī)律與本文中的結果基本一致，尤其是θ=0°的點，兩者得到的σθ/σ值完全相同.但是本文與文獻[11]的結果最大的不同在于，圖3c中最大切向應力的數(shù)值是57.0，發(fā)生在θ=0.5°的點，但是這個最大值的點卻被文獻[11]忽略了.

將這2種各向異性材料與各向同性材料進行對比，如圖2和圖3，發(fā)現(xiàn)正交各向異性板孔邊應力分布的一些特殊規(guī)律：對于各向同性材料，在2種單軸拉伸荷載作用下，孔邊的最大應力集中發(fā)生在θ=0°的尖點位置A，當σx∞=σ單獨作用時，在A點產生的是最大壓應力，而當σy∞=σ單獨作用時，在A點產生的是最大拉應力；但是對于正交各向異性材料，孔邊切向應力的最大值可能發(fā)生在尖點位置，如圖3b中的尖點A，也可能發(fā)生在尖點A的臨近位置，如圖3c中θ=0.5°的點.

除此之外，對于正交各向異性材料，當荷載沿x軸方向作用時(σx∞=σ)，在與x軸相交的孔邊2個點(A點和C點)的切向應力均為壓應力，且其大小都為-σ，如圖2b、2c；當荷載沿y軸方向作用時(σy∞=σ)，在與y軸相交的孔邊2個點的切向應力也均為壓應力，且大小也都是-σ，如圖3b、3c.而對于各向同性材料，在2種荷載作用下都沒有這樣的規(guī)律.但是3種材料的相同點是，當荷載作用方向與孔邊尖點A的指向(x軸方向)垂直時，如圖3，最大的切向應力集中都是發(fā)生在孔邊的A處尖點或其臨近的孔邊點.

2.2 不規(guī)則孔形的孔外域切向應力分析

觀察圖3發(fā)現(xiàn)，在沿y軸方向荷載(σy∞=σ)的作用下，不規(guī)則孔形在θ=0°的孔邊點其切向應力集中非常明顯，為了進一步研究孔外域的應力分布情況，圖4中給出了在θ=0°、180.0°的板內，3種材料結構所對應的切向應力與外荷載的比值σθ/σ，這在文獻[11]中并沒有研究.

由圖4可以看出：各向同性材料最大切向應力產生在孔邊(圖4a，ρ=1.000 0)，而對于正交各向異性材料則不一定，纖維角度按[0°/-90.0°]s布置時，最大切向應力產生在孔邊(圖4b，ρ=1.000 0)，纖維角度按[45.0°/-45.0°]s布置時，最大切向應力產生在非常臨近孔邊的外部區(qū)域，如圖4c中可以清晰地看出：在ρ=1.002 4對應的孔邊右側點，σθ/σ達到最大值34.1，而在ρ=1.000 0的孔邊點，σθ/σ只有29.2；在ρ=1.203 0對應的孔邊左側點，σθ/σ達到最大值2.1，而在ρ=1.000 0的孔邊點，σθ/σ為2.0.

a 各向同性材料

b G-E[0°/-90.0°]s

c G-E[45.0°/-45.0°]s圖4 沿y軸單向拉伸荷載作用下(σy∞=σ)不規(guī)則孔形在θ=0°、180.0°板內點的切向應力分布

Fig.4Distributionofnormalizedtangentialstressesinθ=0°and180.0°forirregularshapedholewithuniaxialloadingσy∞=σ

2.3 正六角形孔的孔邊切向應力分析

所研究的正六角形孔的映射函數(shù)系數(shù)Ck為[13]：C5=0.066 7，C11=0.010 1，C17=0.003 6，C23=0.001 8，C29=0.001 0，其他Ck均為零.材料參數(shù)的取值與上文中不規(guī)則孔完全相同，故特征方程的根也相同.R的值仍為1.0.將各向同性材料與上述2種正交各向異性材料進行對比，圖5和圖6分別用虛線繪出了沿x軸或y軸荷載作用下(σx∞=σ，σy∞=σ)正六角形孔的孔邊切向應力與外荷載的比值σθ/σ的分布.

分析圖5和圖6可以看出：各向同性材料的孔邊最大切向應力發(fā)生在尖點處，如圖5a中θ=60.0°的B點和圖6a中θ=0°的A點.但是對于正交各向異性材料，孔邊產生最大切向應力的位置既可能是在尖點處，如圖6b的尖點A處，也可能是在尖點的附近點，如圖5b中尖點B附近θ=62.0°的點、圖5c中尖點B附近θ=59.5°的點、圖6c中尖點A附近θ=2.5°的點.除此之外，分析圖6還可以看出，當荷載作用方向與孔邊尖點A的指向(x軸方向)垂直時，3種材料最大切向應力都發(fā)生在孔邊A處尖點或者其臨近的孔邊點.這與圖3中不規(guī)則孔形的規(guī)律完全相同.

將圖2和圖3中不規(guī)則孔形的孔邊應力分布與圖5、圖6中正六角形孔的結果進行對比分析發(fā)現(xiàn)：當3種材料結構的平板受沿y軸方向的單向拉伸時，在y軸與孔邊界的交點及其附近的孔邊區(qū)域，切向應力均為壓應力，尤其對于2種正交各向異性板，在y軸與孔邊界的交點，即θ=90.0°的點，切向應力值都是-σ；同樣，當荷載是沿x軸方向拉伸時，在x軸與孔邊界的交點及其附近的孔邊區(qū)域，切向應力也均為壓應力，而2種正交各向異性板中x軸與孔邊界的交點(θ=0°的點)切向應力值也都是-σ.

2.4 孔邊最大切向應力點位置的確定

為了進一步研究正交各向異性板在y方向荷載σy∞=σ單獨作用下孔邊最大切向應力發(fā)生的位置，針對上述2種孔形，計算其他纖維角度工況下孔邊的切向應力與外荷載的比值σθ/σ.不同纖維角度條件下，孔邊產生最大切向應力的位置(θ角)及其最大應力值(σθ/σ)max與尖點(θ=0°)處應力值的對比見表1.

a 各向同性材料

b G-E[0°/-90.0°]s

c G-E[45.0°/-45.0°]s圖5 帶有正六角形孔的3種材料結構在σx∞=σ單獨作用下的孔邊切向應力分布

Fig.5Distributionofnormalizedtangentialstressesforthreekindsofplateswithhexagonholeunderuniaxialloadingσx∞=σ

a 各向同性材料

b G-E[0°/-90.0°]s

c G-E[45.0°/-45.0°]s圖6 帶有正六角形孔的3種材料結構在σy∞=σ單獨作用下的孔邊切向應力分布

Fig.6Distributionofnormalizedtangentialstressesforthreekindsofplateswithhexagonholeunderuniaxialloadingσy∞=σ

表1纖維角度不同時整體坐標系下的復根μ1，μ2以及孔邊最大切向應力(σθ/σ)max與尖點處(θ=0°)應力的對比

Tab.1Thecomplexrootsμ1,μ2intheglobalcoordinatesandthecomparisonofthemaximumnormalizedtangentialstressσθ/σwiththestressinthesharpcorner(θ=0°)fordifferentfiberorientationangles

纖維角度μ1μ2不規(guī)則孔正六角形孔0°時(σθ/σ)(σθ/σ)maxθ/(°)0°時σθ/σ(σθ/σ)maxθ/(°)[0°/-90.0°]s3.6404i0.2747i108.9108.9015.615.60[5.0°/-85.0°]s-0.9732+3.3304i0.0808+0.2766i100.5108.10.0414.415.40.27[10.0°/-80.0°]s-1.5300+2.6583i0.1626+0.2826i82.1105.80.0911.915.00.55[15.0°/-75.0°]s-1.6823+1.9994i0.2464+0.2928i64.2102.00.139.514.20.83[20.0°/-70.0°]s-1.6183+1.4961i0.3332+0.3080i50.796.90.187.713.31.14[25.0°/-65.0°]s-1.4719+1.1418i0.4242+0.3290i41.690.60.236.512.01.48[30.0°/-60.0°]s-1.3058+0.8960i0.5207+0.3573i35.683.30.285.710.61.84[35.0°/-55.0°]s-1.1443+0.7236i0.6243+0.3948i31.875.20.345.29.12.20[40.0°/-50.0°]s-0.9952+0.6005i0.7367+0.4445i29.866.50.404.97.52.47[45.0°/-45.0°]s-0.8597+0.5109i0.8597+0.5109i29.257.00.504.86.02.50

圖7給出了纖維的旋轉方向和其中一種工況([10.0°/-80.0°]s)的示意圖.由于孔形和荷載的對稱性，纖維在[0°/-90.0°]s和[45.0°/-45.0°]s范圍內的角度變化包含了纖維發(fā)生旋轉時所有可能的布置情況，而不同纖維角度的布置對應的整體坐標系下特征方程的根μ1、μ2可以根據式(14)求出，結果見表1.

由表1可見：對于2種孔形，僅當y方向作用有荷載時，纖維角度按[0°/90.0°]s布置，最大的切向應力產生在尖點A(θ=0°)；但當纖維角度稍微發(fā)生旋轉時，最大切向應力的發(fā)生位置逐漸偏離A點，其中，不規(guī)則形狀孔的偏離幅度較小，最遠達到θ=0.5°的點，而正六角形孔的偏離幅度最遠達到θ=2.5°的點；除此之外，隨著纖維角度從[0°/90.0°]s向[45.0°/-45.0°]s旋轉，2種孔形對應的最大切向應力值逐漸減小.無論是不規(guī)則形狀孔還是正六角形孔，正交各向異性板纖維角度的改變對孔邊以及孔外域的切向應力集中都會產生很大的影響，在以后的研究中，可以針對纖維角度進行優(yōu)化設計的研究，以使孔邊和孔外域都獲得較小的應力集中.

圖7 纖維角度按[10.0°/-80.0°]s布置的示意圖

Fig.7Sketchofthefiberorientationanglesat[10.0°/-80.0°]s

3 結論

(1) 在單向均布荷載作用下，各向同性材料孔邊產生最大切向應力的位置是尖點處，而正交各向異性材料則不一定.

(2) 當單向均布荷載的作用方向與孔邊尖點的指向垂直時，纖維角度按[0°/-90.0°]s布置會使最大的切向應力產生在尖點處(θ=0°)；但是隨著纖維角度的旋轉，最大應力點的位置逐漸偏離尖點，最大切向應力的數(shù)值也逐漸減小.

(3) 各向同性材料產生最大切向應力的位置在孔邊，而正交各向異性材料則不一定，當單向拉伸荷載的方向與尖點指向垂直時，最大切向應力可能是在孔邊(纖維按[0°/-90.0°]s角度布置)，也可能是在孔外部的臨近區(qū)域(纖維按[45.0°/-45.0°]s角度布置).

(4) 正交各向異性板切向應力的集中程度也與材料的纖維方向有很大的關系，通過調整材料纖維的方向可以獲得孔邊或者孔外區(qū)域較小的應力集中.

(5) 當無限平板中均布的單向拉伸荷載σ沿某坐標軸方向作用時，該坐標軸與孔邊界的交點及其附近的孔邊區(qū)域，切向應力均為壓應力.尤其對于正交各向異性板，在該坐標軸與孔邊界的交點，切向應力均為-σ.

參考文獻:

[1] LEKHNITSKII S G. Anisotropic plates[M]. New York: Gordon and Breach, 1968.

[2] LEKHNITSKII S G. Theory of elasticity of an anisotropic body[M]. Moscow: Mir Publishers, 1981.

[3] DAOUST J, HOA S V. An analytical solution for anisotropic plates containing triangular holes[J]. Composite Structures, 1991, 19(2): 107.

[4] UKADGAONKER V G, RAO D K N. Stress distribution around triangular holes in anisotropic plates[J]. Composite Structures, 1999, 45(3):171.

[5] SHARMA D S. Stress distribution around circular/elliptical/ triangular holes in infinite composite plate[J]. Engineering Letters, 2012, 20(1): EL_20_1_01.

[6] REZAEEPAZHAND J, JAFARI M. Stress analysis of composite plates with non-circular cutout[J]. Key Engineering Materials, 2008, 385-387:365.

[7] RAO D K N, BABU M R, REDDY K R N,etal. Stress around square and rectangular cutouts in symmetric laminates[J]. Composite Structures, 2010, 92(12):2845.

[8] SHARMA D S. Stresses around polygonal hole in an infinite laminated composite plate[J]. European Journal of Mechanics-A/Solids, 2015, 54:44.

[9] 李成, 鄭艷萍. 不同孔口形狀對含孔復合材料板孔邊應力狀態(tài)影響的研究[J]. 工程力學, 2007, 24(10):19.

LI Cheng, ZHENG Yanping. Influence of different orifice figure on hole-edge stress of composite material plate with hole [J]. Engineering Mechanics, 2007, 24(10): 19.

[10] PATEL N P, SHARMA D S. On the stress concentration around a polygonal cut-out of complex geometry in an infinite orthotropic plate[J]. Composite Structures, 2017, 179: 415.

[11] UKADGAONKER V G, KAKHANDKI V. Stress analysis for an orthotropic plate with an irregular shaped hole for different in-plane loading conditions—Part 1[J]. Composite Structures, 2005, 70(3):255.

[12] LU A Z, ZHANG N, ZHANG X L,etal. Analytic method of stress analysis for an orthotropic rock mass with an arbitrary-shaped tunnel[J]. International Journal of Geomechanics, 2015, 15(4): 04014068.

[13] 呂愛鐘, 張路青. 地下隧洞力學分析的復變函數(shù)方法[M]. 北京：科學出版社, 2007.

Lü Aizhong, ZHANG Luqing. Complex function method on mechanical analysis of underground tunnel[M]. Beijing: Science Press, 2007.

[14] 薩文. 孔附近的應力集中[M]. 盧鼎霍,譯. 北京：科學出版社, 1958.

SAVIN G N. Stress concentration around holes[M]. Translated by LU Dinghuo. Beijing: Science Press, 1958.

[15] 陳子蔭. 圍巖力學分析中的解析方法[M]. 北京：煤炭工業(yè)出版社, 1994.

CHEN Ziyin. Analytical method of rock mechanics analysis[M]. Beijing: China Coal Industry Publishing House, 1994.