張營川 馬 駿

大連理工大學(xué)船舶工程學(xué)院，遼寧大連 116024

1 引 言

螺旋槳振動噪聲，尤其是共振引起的異常噪聲——“唱音”，是艦船三大噪聲源之一。以水下隱蔽性為主要特征的潛艇和聲自導(dǎo)魚雷，為避免過早被敵方發(fā)現(xiàn)，又能及早發(fā)現(xiàn)敵方，必須具有很低的自噪聲，這是提高潛艇、聲自導(dǎo)魚雷隱蔽性和戰(zhàn)術(shù)性能的重要措施。因此，研究螺旋槳“唱音”具有重要意義。

螺旋槳“唱音”是由湍流激勵引起螺旋槳槳葉結(jié)構(gòu)共振而引發(fā)的非線性聲學(xué)現(xiàn)象，這一現(xiàn)象的本質(zhì)可以歸結(jié)到螺旋槳的旋渦發(fā)放問題中，也就是說，在螺旋槳的旋轉(zhuǎn)過程中，槳葉的隨邊產(chǎn)生了旋渦，當(dāng)旋渦發(fā)放的頻率與槳葉的固有頻率接近時，螺旋槳槳葉由于共振便開始發(fā)出清脆的鳴音。雖然對于螺旋槳噪聲的研究已開展很多［1-2］，但至目前為止，對于螺旋槳“唱音”現(xiàn)象的研究還只停留在依靠模型實驗與經(jīng)驗公式階段［3-5］，而且只能通過分析測量得到的螺旋槳的聲學(xué)數(shù)據(jù)來研究避免“唱音”現(xiàn)象的發(fā)生。

近年來，螺旋槳推進領(lǐng)域的數(shù)值研究有所發(fā)展，出現(xiàn)了較多文獻［6-9］，但對于螺旋槳“唱音”問題的數(shù)值研究還少有涉及。如果能夠從計算機輔助工程的角度出發(fā)，建立大小合理的計算模型并能得到工程上精確的螺旋槳槳葉在不同來流速度和不同轉(zhuǎn)速下的旋渦發(fā)放頻率，就可以在螺旋槳的設(shè)計階段通過輸入螺旋槳的模型參數(shù)及轉(zhuǎn)速和航速來計算出螺旋槳是否會發(fā)生“唱音”現(xiàn)象，這是非常具有實際應(yīng)用意義的，其結(jié)果不僅可以減少模型實驗的時間和費用，而且還可避免在船舶下水之后因發(fā)生“唱音”現(xiàn)象而再來修改槳葉隨邊形狀等一系列耗費人力物力的工作。但在當(dāng)前的計算能力下，不做任何簡化，詳細(xì)地建立螺旋槳與周圍流體模型進行流固耦合計算還存在困難。

因此，本文從數(shù)值模擬的角度在螺旋槳 “唱音”問題上做了一些初步研究，即基于LS-DYNA軟件以及任意拉格朗日—歐拉（ALE）方法對螺旋槳和周圍流體建立了有限元模型，通過適當(dāng)?shù)暮喕?，觀察到了槳葉后方旋渦的運動特性，并給出了對應(yīng)的解釋，最后計算并分析了槳葉的結(jié)構(gòu)響應(yīng)及相關(guān)影響因素。這些結(jié)論對于后續(xù)開展螺旋槳渦激共振和“唱音”現(xiàn)象方面的數(shù)值研究具有一定的參考價值。

2 ALE方法的基本公式

ALE方法是描述流體運動的兩種基本方法。在拉格朗日描述下，網(wǎng)格隨流體一起運動，網(wǎng)格點的速度與當(dāng)?shù)亓黧w微團的速度相同。在歐拉描述下，網(wǎng)格的空間位置固定，即網(wǎng)格點的速度為零，流體微團穿過網(wǎng)格單元所構(gòu)成的控制體。文獻［10］的方法是將這兩種基本方法統(tǒng)一起來，允許網(wǎng)格以任意速度運動。當(dāng)網(wǎng)格速度為零時，為歐拉方法，當(dāng)網(wǎng)格速度等于流體運動速度時，則為拉格朗日方法。

在ALE方法的描述中，引入了拉格朗日和歐拉坐標(biāo)之外的第3個任意參照坐標(biāo)。與參照坐標(biāo)相關(guān)的材料微商可采用下式描述：

式中，f為ALE方法中與坐標(biāo)和時間相關(guān)的任意變量；Xi為拉格朗日坐標(biāo)；xi為歐拉坐標(biāo)；wi為相對速度。相對速度通過下式表達(dá)：

式中，v為物質(zhì)速度；u為網(wǎng)格速度。

因此，用材料時間導(dǎo)數(shù)和參照幾何構(gòu)形的時間導(dǎo)數(shù)兩者之間的替換關(guān)系可以推導(dǎo)出所需的ALE方程。

ALE算法的控制方程可由下列守恒方程給出：

1）質(zhì)量守恒方程

式中，ρ為流體密度。

2）動量守恒方程

控制固定域上的牛頓流體流動問題的增強形式由控制方程和對應(yīng)的初始及邊界條件組成，控制流體問題的方程是Navier-Stokes方程的ALE描述：

式中，bi為單位質(zhì)量的體積力分量；σij為Cauchy應(yīng)力張量。

Cauchy應(yīng)力張量能夠以Stokes本構(gòu)方程的形式表示，具體表達(dá)式如下：

式中，p為壓力；δij為kronec ker δ函數(shù)；μ為流體的動力粘性系數(shù)。

式（5）可以與下列邊界條件聯(lián)立求解：

式（6）、式（7）中，nj為邊界外法線單位向量；Ui0為指定初始邊界。

并且，對于計算域有下式成立：

式（8）、式（9）中，Γ 為計算域的完整邊界；Γ1，Γ2為計算域的部分邊界。

假設(shè)整個計算域的速度場在t＝0時刻已知，即

式中，E為內(nèi)能。

3 平面結(jié)構(gòu)旋渦發(fā)放現(xiàn)象數(shù)值模擬

1）圓柱剖面

螺旋槳在旋轉(zhuǎn)時其槳葉和流體的相互作用是非常復(fù)雜的，在模型實驗中也難以觀察到槳后流場情況，直接模擬這一流固耦合過程難免出現(xiàn)誤差。因此，首先應(yīng)針對最簡單的平面圓柱在來流下的旋渦發(fā)放現(xiàn)象進行模擬，然后根據(jù)結(jié)果調(diào)整一系列的計算參數(shù)，只有這樣，才能使后續(xù)的螺旋槳槳葉及流體的模擬盡量接近實際情況。

圖1所示為平面圓柱和流體的有限元模型示意圖以及局部網(wǎng)格放大圖。由于采用了ALE方法處理圓柱與流體的流固耦合問題，因此在建模過程中，圓柱和流體可以分別建模和劃分網(wǎng)格，相互之間互不干擾，網(wǎng)格可以重疊。在數(shù)值求解過程中，各時刻交替地求解流體和固體區(qū)域，并且在交替過程中通過罰函數(shù)方程在耦合界面進行有關(guān)物理量的傳遞，從而達(dá)到不同求解域的相互耦合。

圓柱的直徑為B=0.001m，來流速度為v=4.5 m／s，Strouhal數(shù)取為 St=0.2，可以得出該圓柱的旋渦發(fā)放頻率ns約為：

則圓柱的旋渦發(fā)放周期約為0.001 s。

圖2為數(shù)值模擬得到的圓柱在垂直來流方向的受力曲線。流體流經(jīng)圓柱時，圓柱每發(fā)放出一個旋渦，圓柱在垂直來流方向的受力便會改變一次大小和方向，因此，通過圓柱在垂直來流方向的受力曲線，可以計算出圓柱的旋渦發(fā)放頻率約為0.001 s，這個結(jié)果與公式計算得到的結(jié)果基本吻合。圖3給出了圓柱旋渦發(fā)放時，流體單元的應(yīng)變云圖以及整個流體域的速度云圖。圖4給出了圓柱旋渦發(fā)放時，整個流體域的流體速度云圖（含速度方向）以及局部放大圖像。

盡管如此，數(shù)值計算結(jié)果還是存在著不足，這主要是由于在建模過程中，考慮到實際的計算能力，為了減小整個有限元模型的網(wǎng)格數(shù)量，一方面因是采用單層三維網(wǎng)格來研究厚度很薄的圓柱模型與來流的相互作用，因此沒有建立有一定長度的圓柱模型參與模擬；另一方面，本文對流體區(qū)域的網(wǎng)格大小進行了多次嘗試，在網(wǎng)格較大時，得出的旋渦發(fā)放周期和圓柱受力都不準(zhǔn)確，而在網(wǎng)格較小時，能得到相對較好的計算結(jié)果，但網(wǎng)格的細(xì)化使計算時間大幅增加。所以在計算條件允許的前提下，根據(jù)圓柱直徑為0.001m的情況，流體網(wǎng)格大小應(yīng)設(shè)定為0.00005m。

2）機翼剖面

在上述計算的基礎(chǔ)上，本文數(shù)值模擬了機翼剖面的旋渦發(fā)放現(xiàn)象。選取的翼型為NACA 0012，翼展 0.6m，弦長 0.1 m。 計算模型右側(cè)邊界面施加了流入邊界條件，來流速度為4.5 m／s，左側(cè)邊界面施加了流出邊界條件，其他邊界面施加的是無反射邊界條件，整個水域施加了初始流入速度，可以保證流場的連貫；時間步長Δt受到Courant穩(wěn)定條件的限制取得非常小，在計算過程中，還采用了LS-DYNA軟件中的質(zhì)量縮放方法來節(jié)省計算時間。在計算中，考慮了單層機翼剖面受到流體力作用后自身的應(yīng)力和應(yīng)變情況。

圖5給出了機翼剖面開始旋渦發(fā)放以后的流體速度等值線圖（左）和速度矢量圖（右），并且在需要放大的位置進行了標(biāo)記。圖6為速度矢量圖的放大圖，在圖中用圓圈標(biāo)記了旋渦的位置以及旋渦旋轉(zhuǎn)的方向，清晰地展示了單層機翼剖面在來流速度較低時的旋渦發(fā)放現(xiàn)象。

從單層機翼剖面的尾部最末端選取一個節(jié)點，圖7中的曲線表示的就是該節(jié)點垂直來流方向的位移隨時間的變化情況。該條曲線可以用來讀取機翼剖面的旋渦發(fā)放頻率等重要信息。除流體數(shù)據(jù)外，數(shù)值模擬還得到了機翼剖面在流體力作用下自身的應(yīng)力分布情況，從圖8中可以看出，顏色深的部分為機翼剖面的高應(yīng)力區(qū)域，而顏色較深區(qū)域中的結(jié)構(gòu)單元應(yīng)力比較小。

需要說明的是，由于控制旋渦發(fā)放的是機翼尾部最后端的寬度 （與機翼的整體尺寸相比非常?。?，因此流體網(wǎng)格必須按照此寬度進行設(shè)定才能準(zhǔn)確地與固體網(wǎng)格耦合，從而導(dǎo)致整個計算模型的單元數(shù)量非常多。另外，求解含有結(jié)構(gòu)應(yīng)力應(yīng)變的流固耦合問題要比單獨求解剛性固體的流固耦合問題復(fù)雜得多，需要的計算時間也更長。為了保證能在短時間內(nèi)模擬出機翼剖面的旋渦發(fā)放現(xiàn)象，只能將流體網(wǎng)格適當(dāng)放大，但由此會導(dǎo)致旋渦發(fā)放頻率與實際有所差別。要得到比較精確的翼型剖面旋渦發(fā)放頻率，就需要將本文的模型更加細(xì)化。

4 固定螺旋槳槳葉旋渦發(fā)放現(xiàn)象的數(shù)值模擬

本文選定四葉的MAU4-55螺旋槳作為研究對象。雖然應(yīng)用本文方法不僅可以求出槳葉后方的流場特性，而且還可計算出槳葉內(nèi)部應(yīng)力等結(jié)構(gòu)力學(xué)性能，但是數(shù)值模擬螺旋槳整槳的流固耦合現(xiàn)象需要相當(dāng)長的計算時間，為了提高計算速度，本文假設(shè)：1）只研究螺旋槳的一個槳葉；2）將槳葉根部固定，槳葉不發(fā)生運動，在槳葉導(dǎo)邊前方施加平行來流，來流速度為4.5 m／s，從槳葉導(dǎo)邊向隨邊方向流動；3）調(diào)整槳葉的位置，使槳葉的葉面盡量與來流方向平行；4）槳葉與流體網(wǎng)格較大，以便能在短時間內(nèi)模擬出三維的漩渦發(fā)放現(xiàn)象。

圖9給出了簡化后的流固耦合有限元模型，其中深色網(wǎng)格為槳葉固體網(wǎng)格，淺色為流體網(wǎng)格。

通過數(shù)值模擬固定槳葉的旋渦發(fā)放，給出流體單元的壓力云圖如下。

1）槳葉后方的單個旋渦

如圖10所示，本文的數(shù)值模型和計算方法能夠模擬出固定的螺旋槳槳葉在一定來流速度下的旋渦生成現(xiàn)象，而且展示了槳葉后部旋渦的拉扯破裂情況。經(jīng)分析得出，流體在流經(jīng)槳葉葉面達(dá)到槳葉隨邊時，產(chǎn)生了旋渦發(fā)放，生成一個旋渦。由于槳葉根部結(jié)構(gòu)與槳葉頂部結(jié)構(gòu)的偏轉(zhuǎn)方向不同（圖9中右圖），槳葉根部對應(yīng)的旋渦下部與槳葉頂部對應(yīng)的旋渦上部的運動方向相反，整個旋渦在繼續(xù)運動過程中，使得其因自身不同部分的拉扯而最終導(dǎo)致旋渦斷裂。這一點與圓柱形狀結(jié)構(gòu)的旋渦發(fā)放明顯不同。

2）槳葉后方的2個旋渦

隨著時間的推移，槳葉后方不斷產(chǎn)生新的旋渦。如圖11所示，當(dāng)槳葉后方存在2個旋渦時，一方面，這兩個旋渦在運動過程中可能會發(fā)生接觸，從而產(chǎn)生相互作用，兩個旋渦形狀也會發(fā)生一定的變化；另一方面，由于每個旋渦都有各自的流體力學(xué)特性，因此有的旋渦是在上下不同部分的拉扯下發(fā)生破裂，有的旋渦則不會破裂，或者破裂之后由于吸力還能夠再次組合成一個完整的旋渦。

3）槳葉后方的3個旋渦

如圖12左圖所示，從圖中可以清晰地觀察到3個獨立的旋渦，這表明本文的模型已經(jīng)完全可以模擬出固定槳葉后方的旋渦，從開始生成一直到移動出計算模型邊界的整個運動過程，給出了槳葉后方流場的實時變化情況。在右圖中，3個獨立的旋渦在運動過程中發(fā)生了接觸，進而相互之間的邊界不再清晰，而且導(dǎo)致3個旋渦的流體壓力整體下降。這可以從能量守恒的角度加以解釋，即由于在旋渦的相互接觸過程與流場的紊亂過程中都消耗了旋渦的能量，因此從流體的壓力云圖上便顯示出旋渦壓力的下降，即從高壓力變?yōu)榈蛪毫Α?/p>

接下來，將對槳葉在流體作用下的內(nèi)部結(jié)構(gòu)力學(xué)性能進行分析。

1）槳葉葉梢與槳葉根部的響應(yīng)比較

選取節(jié)點進行速度響應(yīng)比較，數(shù)值模擬結(jié)果如圖13所示 （Z方向速度即為垂直來流方向的速度）。其中，曲線A代表槳葉葉梢節(jié)點的數(shù)值模擬結(jié)果，曲線B代表槳葉根部節(jié)點的數(shù)值模擬結(jié)果。

選取單元進行內(nèi)部應(yīng)力響應(yīng)比較，數(shù)值模擬結(jié)果如圖14所示。其中，曲線A代表槳葉葉梢單元的數(shù)值模擬結(jié)果，曲線B代表槳葉根部單元的數(shù)值模擬結(jié)果。

根據(jù)數(shù)值模擬結(jié)果，在本文的算例中，槳葉葉梢節(jié)點的速度響應(yīng)要大于槳葉葉根節(jié)點的速度響應(yīng)，這主要是因為槳葉的葉根部分比槳葉的葉梢部分更靠近約束條件，而且葉根的厚度也比葉梢的厚度大得多，因此，在槳葉旋渦發(fā)放產(chǎn)生的渦激力作用下，葉根部分的節(jié)點難以發(fā)生運動。

對等截面的懸臂梁來說，靠近約束的位置其單元應(yīng)力較大；而對于變截面的懸臂梁來說其單元應(yīng)力與梁的厚度有很大關(guān)系，厚度大的位置其單元應(yīng)力較小。本文的槳葉模型近似于變截面的懸臂梁，但厚度大的位置卻靠近約束，所以，盡管圖10的數(shù)值模擬結(jié)果顯示槳葉葉根部分的單元應(yīng)力比葉梢部分更大，但這只是針對本文的算例而言，對于具體的槳葉，需具體計算才能得到正確的結(jié)論。

空白組大鼠踝關(guān)節(jié)表面光滑，滑膜細(xì)胞無增生，周圍未見炎癥細(xì)胞浸潤。模型組大鼠關(guān)節(jié)間隙變窄，關(guān)節(jié)軟骨退變，局部可見骨質(zhì)破壞，滑膜表面粗糙，滑膜細(xì)胞增生、層次紊亂，滑膜及周圍軟組織可見大量炎癥細(xì)胞浸潤，以淋巴細(xì)胞、巨噬細(xì)胞浸潤為主，可見豐富的血管翳形成，纖維母細(xì)胞和脂肪細(xì)胞增生活躍。與模型組比較，各給藥組大鼠關(guān)節(jié)軟骨退變、骨破壞及滑膜細(xì)胞增生明顯減輕，血管翳形成較少，炎癥細(xì)胞數(shù)量減少，僅見少量纖維母細(xì)胞及脂肪細(xì)胞增生，詳見圖2。

2）槳葉隨邊與槳葉葉面中部的響應(yīng)比較

選取槳葉隨邊和葉面中部的節(jié)點進行速度響應(yīng)比較，數(shù)值模擬結(jié)果如圖15所示。同時，選取槳葉隨邊和葉面中部的單元進行有效應(yīng)力響應(yīng)比較，數(shù)值模擬結(jié)果如圖16所示。其中，曲線A代表槳葉隨邊的數(shù)值模擬結(jié)果，曲線B代表槳葉中部的數(shù)值模擬結(jié)果，選取的節(jié)點和單元至槳軸的距離基本相同。

圖中線條表明，由于選取的兩個節(jié)點至槳葉根部約束的距離基本相同，而槳葉隨邊的厚度又小于槳葉中部的厚度，所以，槳葉隨邊的節(jié)點速度響應(yīng)要比槳葉中部的大。另一方面，由于流體流經(jīng)槳葉葉面時，在槳葉隨邊產(chǎn)生了旋渦，旋渦脫落產(chǎn)生的渦激力大部分作用在槳葉的隨邊，而且槳葉隨邊的厚度要小于槳葉中部的厚度，因此，槳葉隨邊單元的有效應(yīng)力響應(yīng)大于槳葉中部單元的有效應(yīng)力響應(yīng)。

5 結(jié) 論

綜上所述，經(jīng)過一系列的假定和簡化，采用數(shù)值模擬的方法研究螺旋槳固定槳葉的旋渦發(fā)放現(xiàn)象是完全可行的，并且觀察得到了3個主要結(jié)論：

1）螺旋槳槳葉根部與頂部結(jié)構(gòu)的偏向不同，在流體流動到槳葉隨邊產(chǎn)生旋渦發(fā)放時，旋渦上部與下部的運動方向也不相同，因此，部分旋渦會在運動過程中發(fā)生破裂，但旋渦破裂是否發(fā)生、破裂發(fā)生后旋渦是否還會重組，這與旋渦本身的流體特性有關(guān)。

2）獨立的旋渦在逐漸遠(yuǎn)離槳葉的過程中，容易與其他旋渦發(fā)生碰撞而導(dǎo)致旋渦形狀的改變和旋渦邊界的破壞，以至整個流場的紊亂，這些都消耗了旋渦的能量，最終降低了旋渦的壓力。

3）與槳葉葉根部分相比，槳葉葉梢部分的節(jié)點速度較大。但對于不同的槳葉，需具體計算才能比較槳葉葉梢與葉根部分的內(nèi)部應(yīng)力大小。與槳葉中部相比，槳葉隨邊部分的節(jié)點速度較大，單元有效應(yīng)力也較大。槳葉內(nèi)部有效應(yīng)力的主要影響因素在于選取的單元所在位置的槳葉厚度，以及單元至槳葉根部的距離和單元至旋渦脫落位置的距離。

上述研究只是螺旋槳“唱音”問題的基礎(chǔ)性工作，要從數(shù)值模擬的角度完整研究螺旋槳“唱音”的機理和各種影響因素，還需從以下幾個方面進行努力：

首先，螺旋槳“唱音”現(xiàn)象的實質(zhì)是槳葉后方的旋渦發(fā)放頻率接近于槳葉的固有頻率，因此，通過數(shù)值模擬的方法求出旋轉(zhuǎn)中的螺旋槳槳葉后方旋渦的發(fā)放頻率是問題的關(guān)鍵。由于受條件所限，本文僅給出了定性的分析。其原因在于，螺旋槳槳葉的隨邊厚度直接影響著旋渦發(fā)放的頻率，而通常槳葉隨邊的厚度相對于整個螺旋槳的大小來說非常小，而在流固耦合數(shù)值模擬過程中，又必須要以隨邊的厚度來決定流體和固體網(wǎng)格的大小。因此，要計算得出較接近于實際情況的槳葉旋渦發(fā)放頻率，整個數(shù)值模型就需要用大量的單元來進行描述。但完成這項工作離不開大型計算設(shè)備的支持，同時對研究者的模型簡化能力也是一種考驗。

其次，本文采用的簡化模型是將槳葉固定，在槳葉導(dǎo)邊前方施加流速來模擬螺旋槳轉(zhuǎn)動時槳葉和流體的相對運動。這種方法雖然反應(yīng)出了槳葉渦激共振問題的本質(zhì)，但相對于完整的螺旋槳旋轉(zhuǎn)過程來說，還不能體現(xiàn)各個槳葉后方所產(chǎn)生的旋渦之間的相互作用現(xiàn)象。

最后，目前螺旋槳“唱音”問題的相關(guān)研究都只考慮了螺旋槳轉(zhuǎn)速與“唱音”現(xiàn)象之間的關(guān)系，沒有加入螺旋槳前方來流的因素。而實際上螺旋槳在靜水中旋轉(zhuǎn)和在不同來流下旋轉(zhuǎn)，槳葉隨邊后方的旋渦發(fā)放方向與頻率都會發(fā)生變化。也就是說，帶入旋渦發(fā)放頻率計算公式中的速度v應(yīng)該是一種合速度，即螺旋槳旋轉(zhuǎn)引起的槳葉隨邊與靜止流體的相對速度和前方來流速度的合速度。雖然由于螺旋槳槳葉形狀不規(guī)則，難以從理論上計算出合速度的大小和方向，但在研究螺旋槳“唱音”問題時，在數(shù)值模擬過程中，必須考慮前方來流與螺旋槳轉(zhuǎn)速的共同作用，才能接近于真實情況。

［1］SEOl H，JUNG B，SUH J C.Prediction of non-cavitating underwater propeller noise ［J］.Journal of Sound and Vibration，2002，257（1）：131-156.

［2］KEHR Y Z，KAO J H.Numerical prediction of the blade rate noise induced by marine propellers［J］.Journal of Ship Research，2004，48（1）：1-14 .

［3］于大鵬，趙德有.螺旋槳鳴音的混沌動力特性研究［J］.振動與沖擊，2009，28（12）：47-52.

［4］李潔雅，歐禮堅，王志勇.對螺旋槳引起的諧鳴振動的探討與實踐［J］.中國修船，2001，15（6）：26-27.

［5］華漢金.螺旋槳鳴音產(chǎn)生機理及防治方法 ［J］.船舶，2002（2）：20-23.

［6］劉丹，陳鳳馨.CFD在計算船舶螺旋槳敞水性能中的應(yīng)用研究［J］.現(xiàn)代制造工程，2010（4）：18-20.

［7］賴華威，劉月琴，吳家鳴.基于CFD方法的螺旋槳性能計算與分析［J］.船海工程，2009，38（4）：131-135.

［8］胡健，黃勝，王培生.螺旋槳水動力性能的數(shù)值預(yù)報方法［J］.江蘇科技大學(xué)學(xué)報，2008，22（1）：1-6.

［9］蔡榮泉，陳鳳明，馮學(xué)梅.采用Fluent軟件的螺旋槳敞水性能計算分析［J］.船舶力學(xué)，2006，10（5）：41-48.

［10］李裕春，時黨勇，趙遠(yuǎn).ANSYS 10.0／LS-DYNA 基礎(chǔ)理論與工程實踐［M］.北京：中國水利水電出版社，2008.