高 峰，趙馮兵

（重慶交通大學(xué) 土木建筑學(xué)院，重慶 400074）

人工邊界從廣義上可分為靜力人工邊界和動(dòng)力人工邊界。靜力人工邊界由來已久，通常有固定邊界、滾軸邊界等。動(dòng)力人工邊界經(jīng)過幾十年的發(fā)展，已形成全局人工邊界和局部人工邊界兩大類別，并應(yīng)用于各自適應(yīng)的工程計(jì)算中［1］。此外，王后裕、陳上明等［2］將透射邊界應(yīng)用于考慮結(jié)構(gòu)與圍巖動(dòng)力相互作用的地下離壁式框架結(jié)構(gòu)地震分析模型，陳萬祥、郭志昆等［3］將粘彈性支承應(yīng)用于低速?zèng)_擊下的邊界梁反應(yīng)，均取得了較好的效果。

動(dòng)力人工邊界發(fā)展到現(xiàn)在已有透射邊界、粘性邊界、粘彈性邊界等幾種類型。Deeks［4］提出粘彈性人工邊界。劉晶波等人［5，6］發(fā)展了二維的粘彈性人工邊界，以后將其發(fā)展為三維時(shí)域粘彈性人工邊界。盧華喜等人［7］考慮土的動(dòng)力非線性性質(zhì)，基于土體等效線性模型中的動(dòng)模量變化模式，推導(dǎo)了適用于土－結(jié)構(gòu)動(dòng)力相互作用的粘彈性人工邊界。劉晶波等人［8，9］將二維粘彈性邊界發(fā)展成一致粘彈性人工邊界及其對應(yīng)的粘彈性邊界單元，后來又推導(dǎo)了三維一致粘彈性人工邊界及等效粘彈性邊界單元。

地下結(jié)構(gòu)動(dòng)力響應(yīng)分析中一般需要考慮初始應(yīng)力場，可稱之為靜－動(dòng)力分析問題，其一般算法為先對模型進(jìn)行靜力分析，求出靜應(yīng)力場，然后以此為初始條件進(jìn)行動(dòng)力分析。目前對線彈性小變形范圍內(nèi)的靜－動(dòng)力分析問題，普遍做法是采用靜力人工邊界和動(dòng)力人工邊界分別對靜力和動(dòng)力問題進(jìn)行計(jì)算，將計(jì)算結(jié)果進(jìn)行疊加得到完整的結(jié)果［1］。

對于涉及非線性或大變形問題的靜－動(dòng)力分析，靜力分析轉(zhuǎn)換到動(dòng)力分析時(shí)人工邊界常用處理方法主要有以下幾種：① 靜力分析中采用滾軸邊界或固定邊界，動(dòng)力分析采用粘彈性邊界、透射邊界、粘性邊界等人工邊界；② 靜力分析和動(dòng)力分析都采用靜－動(dòng)力統(tǒng)一邊界，如粘彈性靜－動(dòng)力統(tǒng)一人工邊界；③ 靜力分析和動(dòng)力分析都采用滾軸邊界或固定邊界。

對第③種方法，由于固定邊界使波動(dòng)全部反射，已有許多文獻(xiàn)證明其具有放大振動(dòng)效應(yīng)的作用，目前已經(jīng)使用得不多。劉晶波等人基于粘彈性動(dòng)力人工邊界和半無限空間中靜力問題的基本解，建立了對動(dòng)力問題和靜力問題均適用的三維粘彈性靜－動(dòng)力統(tǒng)一人工邊界，從而上述第②種方法得以解決［1］。

然而，在使用人工邊界對地下結(jié)構(gòu)進(jìn)行動(dòng)力分析時(shí)，還存在一些問題。如第（1）種方法，由于對靜－動(dòng)力分析中的人工邊界轉(zhuǎn)換原理認(rèn)識(shí)不清，將使計(jì)算結(jié)果嚴(yán)重偏離預(yù)期值，失去了施加動(dòng)力人工邊界的必要性。在第②種方法中，將粘彈性靜－動(dòng)力統(tǒng)一人工邊界應(yīng)用于地下結(jié)構(gòu)的靜力分析時(shí)，其計(jì)算結(jié)果誤差較大。本文將就此兩個(gè)問題進(jìn)行論證和分析，并闡述合理的地下結(jié)構(gòu)靜－動(dòng)力分析人工邊界轉(zhuǎn)換方法。

1 靜力和動(dòng)力有限元分析原理

1.1 靜力分析原理

地鐵等地下工程中初始應(yīng)力場的確定須先計(jì)算未開挖狀態(tài)下圍巖的自重應(yīng)力場，進(jìn)而根據(jù)施工步驟，采用釋放荷載法，計(jì)算出襯砌結(jié)構(gòu)和圍巖的靜應(yīng)力場［10］。許多地下結(jié)構(gòu)的自重應(yīng)力場模型可以假設(shè)為半無限空間體，根據(jù)經(jīng)典圍巖壓力理論和彈性力學(xué)理論，半無限空間體中距地表面任一深度h處的應(yīng)力狀態(tài)可定義為［11，12］：

式中，σV為豎向應(yīng)力，γ為圍巖重度，σH為橫向應(yīng)力，λ為側(cè)壓力系數(shù)。對淺層圍巖，可假設(shè)其為各向同性介質(zhì)，側(cè)壓力系數(shù)可表示如下［11］：

對有限元計(jì)算中所取的有限區(qū)域，可以據(jù)此確定有限域邊界條件。

1.2 靜力分析邊界條件

根據(jù)上述靜力分析方法，在靜、動(dòng)力共同作用問題的計(jì)算中，在進(jìn)行動(dòng)力分析之前須先確定地下結(jié)構(gòu)模型的靜應(yīng)力場。地下結(jié)構(gòu)開挖前可將大地假設(shè)為半無限空間體，其在重力作用下的靜力計(jì)算，根據(jù)對稱性，模型中任一處的水平位移uH=0。因此，計(jì)算該應(yīng)力場時(shí)，有限區(qū)域模型兩側(cè)可用水平約束即法向約束，底部可用全約束或僅約束豎直方向，頂面即地面應(yīng)為自由邊界。

1.3 動(dòng)力分析原理

在動(dòng)力荷載作用下，有限元體系在t+Δt時(shí)刻的運(yùn)動(dòng)平衡方程為：

式中α、β為常數(shù)，可按兩種不同的振動(dòng)頻率下測得的阻尼比ξ加以確定。

1.4 粘彈性人工邊界

粘彈性人工邊界從用途上可分為動(dòng)力人工邊界和靜－動(dòng)力統(tǒng)一人工邊界；從具體實(shí)現(xiàn)方法上可分為彈簧－阻尼器邊界單元和一致粘彈性邊界單元。

實(shí)際應(yīng)用中模型邊界單元的物理力學(xué)參數(shù)由其相鄰的圍巖介質(zhì)材料決定，具體計(jì)算公式參考相關(guān)文獻(xiàn)［1，5－9］。結(jié)合有限元法具體實(shí)現(xiàn)時(shí)，此處采用一致邊界，即在人工邊界處設(shè)置一層與主計(jì)算區(qū)域單元有著相同形函數(shù)的單元，其等效剪切模量、等效彈性模量和等效阻尼系數(shù)表達(dá)形式如下：

其中等效阻尼系數(shù)取的是法向和切向的平均值。式中h為等效邊界單元厚度，為等效泊松比，E、G、μ和ρ分別為介質(zhì)彈性模量、剪切模量、泊松比和質(zhì)量密度。R為波源至人工邊界點(diǎn)的距離；cs和cp分別為S波和P波波速；αT與αN分別為切向與法向粘彈性人工邊界參數(shù)，具體取值情況見表1。等效泊松比可按如下取值

表1 粘彈性動(dòng)力人工邊界中參數(shù)α的取值Tab.1 The values of α in viscous-spring artificial boundary for dynamic analysis

表2 粘彈性靜－動(dòng)力統(tǒng)一人工邊界中參數(shù)α的取值Tab.2 The values of α in unified viscous-spring artificial boundary for static and dynamic analysis

靜－動(dòng)力統(tǒng)一人工邊界是在上述動(dòng)力人工邊界的基礎(chǔ)上對人工邊界參數(shù)α進(jìn)行調(diào)整，其具體取值見表2。其中參數(shù)α*具體計(jì)算公式：

式中μ為圍巖泊松比，R為荷載作用點(diǎn)到人工邊界點(diǎn)的距離，d為位置坐標(biāo)。對于底面人工邊界，式中d取荷載作用點(diǎn)至邊界單元的水平距離；對于側(cè)面人工邊界，式中d取荷載作用點(diǎn)至邊界單元的垂直距離。

2 兩種地下結(jié)構(gòu)靜－動(dòng)力人工邊界轉(zhuǎn)換方法邊界轉(zhuǎn)換原理分析

動(dòng)力計(jì)算的初始時(shí)刻，未施加動(dòng)力荷載前，整個(gè)模型各點(diǎn)加速度均為零。地下結(jié)構(gòu)靜－動(dòng)力分析中，要求動(dòng)力分析以靜力分析結(jié)果為初始條件，則模型在動(dòng)力計(jì)算初始時(shí)刻處于靜力平衡狀態(tài)。因此人工邊界的轉(zhuǎn)換可從靜力學(xué)方面討論。

在引言中已經(jīng)介紹了地下結(jié)構(gòu)靜－動(dòng)力分析中人工邊界的常用三種處理方法，下面對第①和第②種方法的人工邊界轉(zhuǎn)換原理進(jìn)行分析。

第①種人工邊界轉(zhuǎn)換方法（方法①）：

該方法靜力分析中采用滾軸邊界或固定邊界（可稱之為傳統(tǒng)靜力邊界條件），動(dòng)力分析中采用粘彈性動(dòng)力人工邊界。一般計(jì)算過程如表3所示，具體為：先采用靜力計(jì)算方法按釋放荷載法計(jì)算地下結(jié)構(gòu)的開挖過程得到初始靜應(yīng)力場等結(jié)果；然后以此靜力計(jì)算結(jié)果為初始條件進(jìn)行動(dòng)力計(jì)算。

表3 方法①靜動(dòng)力邊界條件和荷載條件表Tab.3 Boundaries and loads of the static analysis and the dynamic analysis in method①

下面討論動(dòng)力計(jì)算初始條件中所應(yīng)輸入的靜力計(jì)算結(jié)果。

靜力計(jì)算中總體剛度平衡方程可表為：

式中［K］為剛度矩陣；{F}為等效節(jié)點(diǎn)力向量；{FA}為靜力計(jì)算中所施加的全部外荷載引起的等效節(jié)點(diǎn)力；{FR}為全部約束反力向量，須經(jīng)此方程引進(jìn)邊界條件求出{u}后，才能求得。

靜力計(jì)算中可同時(shí)算得模型靜應(yīng)力場：

式中［D］、［B］和{u}e分別為彈性矩陣、幾何矩陣和單元節(jié)點(diǎn)位移向量。

從此應(yīng)力場的求解過程可以看出，當(dāng)以此應(yīng)力場為初始條件進(jìn)行新的計(jì)算時(shí)，對模型施加引起此應(yīng)力場的荷載和約束是模型保持平衡的必要條件。

由于邊界節(jié)點(diǎn)在約束或指定位移方向不會(huì)對邊界以內(nèi)節(jié)點(diǎn)產(chǎn)生等效節(jié)點(diǎn)力，（11）式可以表示成：

其中：c表示邊界以內(nèi)自由度，b表示邊界自由度，A表示與施加的外荷載相關(guān)，R表示與反力相關(guān)。該式經(jīng)引入邊界條件可得縮減的剛度平衡方程：

由式（14）和式（15）可以解得約束反力引起的節(jié)點(diǎn)荷載向量即約束反力向量為：

去掉約束條件后，將式（16）所示的約束反力代入式（14），仍可解得{ub}為零，且{uc}的解也不變。

這說明對模型施加約束條件與在邊界施加相應(yīng)的約束反力在靜力計(jì)算中是等效的，這使得靜－動(dòng)力分析中以靜力計(jì)算結(jié)果為初始條件時(shí)，可以將約束條件去掉施加以相應(yīng)的約束反力，而不改變模型的靜力平衡狀態(tài)。

綜上可得，在地下結(jié)構(gòu)靜－動(dòng)力分析中，以靜力分析結(jié)果中的應(yīng)力場為初始條件進(jìn)行動(dòng)力分析時(shí)，為保證模型在動(dòng)力計(jì)算初始時(shí)刻為靜力平衡狀態(tài)，可在輸入靜應(yīng)力場的同時(shí)輸入引起此應(yīng)力場的荷載和約束；轉(zhuǎn)換邊界時(shí)，須去掉約束條件，只要同時(shí)施加以相應(yīng)的約束反力即可。

方法①在應(yīng)用中經(jīng)常因?yàn)槲磳?dòng)力計(jì)算初始條件完整輸入模型而導(dǎo)致不合理的結(jié)果，未輸入靜力計(jì)算約束反力就是典型的一種。圖1就是在這種計(jì)算條件下對某地下結(jié)構(gòu)彈性模型在兩種荷載條件下進(jìn)行靜－動(dòng)力分析所得的模型中某一點(diǎn)的動(dòng)力反應(yīng)。圖中振動(dòng)荷載是指如圖8所示的列車振動(dòng)荷載，零動(dòng)力荷載是指動(dòng)力分析中所施加的振動(dòng)荷載（或動(dòng)力荷載）為零。

從圖1（a）和圖1（b）可以看出，該點(diǎn)的位移和應(yīng)力均呈現(xiàn)首先大幅波動(dòng)，然后逐漸趨于穩(wěn)定的情形。從圖1（c）和圖1（d）可以看出，該點(diǎn)的位移和應(yīng)力在沒有任何動(dòng)力荷載情況下仍出現(xiàn)類似波動(dòng)情形。

圖1 某地下結(jié)構(gòu)模型動(dòng)力響應(yīng)Fig.1 Dynamic response in a underground structure modal

模型在任何動(dòng)荷載情況下仍產(chǎn)生較大位移和應(yīng)力，說明模型在零時(shí)刻就受到一個(gè)或一組不平衡力。由此可知模型靜－動(dòng)力轉(zhuǎn)換后未施加靜力分析約束反力將導(dǎo)致不合理的結(jié)果，亦可得此方法中完整施加動(dòng)力計(jì)算初始條件的必要性。

表4 方法②靜動(dòng)力邊界條件和荷載條件表Tab.4 Boundaries and loads of the static analysis and the dynamic analysis in method②

第②種人工邊界轉(zhuǎn)換方法（方法②）：

該方法在靜力計(jì)算與動(dòng)力計(jì)算中均采用粘彈性靜－動(dòng)力統(tǒng)一人工邊界。計(jì)算過程與方法①類似，仍先采用靜力計(jì)算方法按釋放荷載法計(jì)算隧道的開挖過程得到隧道靜應(yīng)力場；然后以此靜應(yīng)力場為初始條件進(jìn)行動(dòng)力計(jì)算。靜力分析和動(dòng)力分析中的荷載和邊界條件如表4所示。

由上述對方法①的分析可知，方法②在計(jì)算過程中邊界條件并未改變，只須按方法①中輸入初始應(yīng)力場的方法輸入初始條件即可。

由此可知，方法②的計(jì)算過程在理論是合理的，但其人工邊界的在具體工程實(shí)例中的靜力計(jì)算精度有待檢驗(yàn)。

3 粘彈性靜－動(dòng)力統(tǒng)一人工邊界的靜力分析精度

上述第②種方法從靜力分析轉(zhuǎn)換到動(dòng)力分析時(shí)靜應(yīng)力場和荷載等的過渡是合理的，下面檢驗(yàn)其相應(yīng)人工邊界的靜力計(jì)算精度。

由于半無限空間體具有經(jīng)典彈性力學(xué)解，即使在巖土力學(xué)中也有理論解，前面已有論述；因此，對半無限空間體的有限元模型進(jìn)行靜力分析，將結(jié)果與理論解進(jìn)行對比，可以檢驗(yàn)其計(jì)算方法的有效性。該模型的理論解可以通過式（1）、式（2）和式（3）算得。

表5 三種工況靜力計(jì)算條件表Tab.5 Static calculation conditions of three conditions

圖2 半無限空間體理論解模型（工況一）Fig.2 Half space modal of the theoretical solution（Condition 1）

圖3 傳統(tǒng)靜力邊界模型（工況二）Fig.3 Static modal of the traditional boundaries（Condition 2）

圖4 粘彈性靜－動(dòng)力統(tǒng)一人工邊界模型（工況三）Fig.4 Static modal of the unified viscous-spring boundary for static and dynamic analysis（Condition 3）

為了比較，進(jìn)行了三種工況的計(jì)算，具體的計(jì)算步驟和計(jì)算條件見表 5。其中工況一為理論解，圍巖性質(zhì)同表6，其計(jì)算模型分別如圖2、圖3和圖4所示。工況三人工邊界彈性模量的計(jì)算采用了動(dòng)力荷載施加點(diǎn)處于模型幾何中心點(diǎn)的假設(shè)。

在模型對稱軸上設(shè)置觀測點(diǎn)，由地面向下每10m設(shè)置一個(gè)觀測點(diǎn)，并將三種工況計(jì)算結(jié)果繪于同一圖中進(jìn)行比較，比較結(jié)果見圖5和圖6。

從圖中可以看出，該有限元模型在施加粘彈性靜－動(dòng)力統(tǒng)一人工邊界時(shí)，其靜力計(jì)算結(jié)果與理論解相比有較大誤差，特別是橫向應(yīng)力結(jié)果誤差較大；而其施加前面所述傳統(tǒng)靜力邊界條件時(shí)，豎向應(yīng)力和橫向應(yīng)力結(jié)果都幾乎與理論解重合。

這說明粘彈性靜－動(dòng)力統(tǒng)一人工邊界應(yīng)用于受重力作用的半無限空間體的有限元模型靜力計(jì)算時(shí)，其誤差偏大。而傳統(tǒng)的靜力邊界即滾軸邊界或固定邊界應(yīng)用于該有限元模型時(shí)精度較好。

4 推薦的地下結(jié)構(gòu)靜－動(dòng)力分析人工邊界轉(zhuǎn)換方法

方法（2）在靜－動(dòng)力分析時(shí)采取了統(tǒng)一的人工邊界，應(yīng)用比較方便，但其相應(yīng)人工邊界靜力分析誤差偏大，從而造成動(dòng)力分析結(jié)果不合理。

由上述對方法（1）的分析得出，較為合理的地下結(jié)構(gòu)的靜－動(dòng)力分析可采用如下計(jì)算步驟：

（1）采用靜力約束條件對模型進(jìn)行靜力分析，計(jì)算結(jié)果應(yīng)包括應(yīng)力場和約束反力。

（2）將靜力模型中的約束條件去掉，施加粘彈性人工邊界。

（3）在原靜力約束被去掉的地方施加約束反力，對模型施加靜力場和靜力荷載（如重力加速度或其他恒力）。

（4）對模型施加動(dòng)力荷載，進(jìn)行動(dòng)力計(jì)算。

模型在動(dòng)力計(jì)算初始時(shí)刻處于靜力平衡狀態(tài)，這是檢驗(yàn)邊界轉(zhuǎn)換合理與否的標(biāo)準(zhǔn)。具體實(shí)現(xiàn)時(shí)可在動(dòng)力計(jì)算中施加零動(dòng)力荷載，并查看模型應(yīng)力位移等的時(shí)程反應(yīng)。

表6 圍巖物理力學(xué)參數(shù)表Tab.6 physical and mechanical parameters of the rock

某雙線鐵路隧道，其凈空跨度 10.2m，為三心圓斷面；襯砌初期支護(hù)15cm，二次襯砌30cm；隧道埋深 45m，Ⅳ類圍巖。采用ANSYS大型有限元分析軟件對該隧道進(jìn)行建模和計(jì)算，模型區(qū)域豎向取100m，橫向取100m。隧道及圍巖物理力學(xué)參數(shù)見表6。

圖7 推薦的方法動(dòng)力分析模型Fig.7 Modal of tunnel for the dynamic analysis in the recommended method

對該隧道進(jìn)行列車荷載作用下的振動(dòng)分析。該算例中的動(dòng)力計(jì)算模型可表示為圖7所示模型，其中FR為靜力分析的約束反力，P（t）為動(dòng)力荷載。

圖8 振動(dòng)荷載時(shí)程曲線Fig.8 Time history curve of the vibration load

計(jì)算分為在動(dòng)力分析中施加振動(dòng)荷載（圖8所示）和零動(dòng)力荷載兩種情況。

圖9（a）和圖9（b）是施加振動(dòng)荷載情況下拱頂?shù)呢Q向位移和第一主應(yīng)力時(shí)程曲線。圖9（c）和圖9（d）是施加零動(dòng)力荷載情況下拱頂?shù)呢Q向位移和第一主應(yīng)力時(shí)程曲線圖。

從圖9（a）和9（b）可以看到，拱頂?shù)呢Q向位移和第一主應(yīng)力均呈現(xiàn)初期波動(dòng)較大，然后逐漸演變?yōu)榉€(wěn)態(tài)振動(dòng)；豎向位移最大振幅約0.045mm，第一主應(yīng)力最大振幅約140Pa。而圖9（c）和圖9（d）顯示在零動(dòng)力荷載作用下，隧道拱頂發(fā)生極小的振動(dòng)，豎向位移最大振幅約1.4×10－6mm，第一主應(yīng)力無振動(dòng)，此可視為數(shù)值計(jì)算誤差，可以忽略不計(jì)。據(jù)此可以推斷動(dòng)力分析初始時(shí)刻模型處于很好的平衡狀態(tài)，從而說明其人工邊界轉(zhuǎn)換更為合理。

圖9 拱頂動(dòng)力響應(yīng)Fig.9 Calculated dynamic response at the vault crown

5 結(jié)論

本文通過將粘彈性人工邊界應(yīng)用于同時(shí)考慮靜力效應(yīng)和動(dòng)力效應(yīng)的地下結(jié)構(gòu)靜－動(dòng)力分析，推薦了一種較為合理的地下結(jié)構(gòu)靜－動(dòng)力分析的人工邊界轉(zhuǎn)換方法，并得到以下結(jié)論：

（1）靜－動(dòng)力分析人工邊界轉(zhuǎn)換時(shí)須保證模型動(dòng)力分析在零時(shí)刻為靜力平衡體；靜－動(dòng)力人工邊界轉(zhuǎn)換方法準(zhǔn)確性，可以通過在動(dòng)力分析中施加零動(dòng)力荷載來進(jìn)行檢驗(yàn)。

（2）在分析過程中均采用粘彈性靜－動(dòng)力統(tǒng)一人工邊界時(shí)，該靜－動(dòng)力分析從靜力狀態(tài)到動(dòng)力分析的過渡是合理的。但該邊界在半無限體有限元模型的靜力計(jì)算中誤差較大，因此，該人工邊界在地下結(jié)構(gòu)靜－動(dòng)力分析中應(yīng)謹(jǐn)慎使用。

［1］劉晶波，李 彬.三維黏彈性靜－動(dòng)力統(tǒng)一人工邊界［J］.中國科學(xué)（E輯），2005，35（9）：966－980

［2］王后裕，陳上明，言志信，等.考慮動(dòng)力相互作用的地下平面鋼框架地震反應(yīng)分析［J］.振動(dòng)與沖擊，2005，24（1）：99－103.

［3］陳萬祥，郭志昆.粘彈性邊界梁在低速?zèng)_擊下的動(dòng)力響應(yīng)分析［J］.振動(dòng)與沖擊，2008，27（12）：69－72，78，108.

［4］ Deeks A J.Randolph M F.Axisymmetric time-domain transmitting boundaries［J］. Journal of Engineering Mechanics，1994，120（1）：25 － 42.

［5］劉晶波，呂彥東.結(jié)構(gòu)－地基動(dòng)力相互作用問題分析的一種直接方法［J］.土木工程學(xué)報(bào)，1998，31（3）：55－64.

［6］劉晶波，王振宇，杜修力，等.波動(dòng)問題中的三維時(shí)域粘彈性人工邊界［J］.工程力學(xué)，2005，22（6）：46－51.

［7］盧華喜，梁平英，尚守平.地基非線性波動(dòng)問題中黏－彈性人工邊界研究［J］.巖土力學(xué)，2005，29（7）：1911－1916.

［8］劉晶波，谷 音，杜義欣.一致粘彈性人工邊界及粘彈性邊界單元［J］.巖土工程學(xué)報(bào)，2006，28（9）：1070－1075.

［9］谷 音，劉晶波，杜義欣.三維一致粘彈性人工邊界及等效粘彈性邊界單元［J］.工程力學(xué)，2007，24（12）：31－37.

［10］高 峰，關(guān)寶樹.深圳地鐵地震反應(yīng)分析［J］.西南交通大學(xué)學(xué)報(bào)，2001，36（4）：355－359.

［11］關(guān)寶樹.隧道力學(xué)概論［M］.成都：西南交通大學(xué)出版社，1993.

［12］徐芝綸.彈性力學(xué)（上冊）［M］.北京：高等教育出版社，2006.

［13］潘昌實(shí).隧道力學(xué)數(shù)值方法［M］.北京：中國鐵道出版社，1995.