時間：2023-06-07 09:09:10
開篇：寫作不僅是一種記錄，更是一種創造，它讓我們能夠捕捉那些稍縱即逝的靈感，將它們永久地定格在紙上。下面是小編精心整理的12篇數值分析，希望這些內容能成為您創作過程中的良師益友，陪伴您不斷探索和進步。
在計算機上運用數值分析解決實際問題的過程為：實際問題→數學模型→數值計算方法→程序設計→上級計算求出結果[1]。數值實驗的設計應當充分體現這一過程，同時也應當充分體現數學建模思想. 而目前地方高校在數值分析課程的教學與實驗中，普遍存在重理論、輕實踐、純數學十足的問題，部分即使避免這些問題，做到利用計算機進行可視化教學和算法編程實踐，但學生也缺乏分析問題、解決問題的能力，因為學生的算法與程序多是百度而來，缺乏思考. 因此，數值分析課程的數值實驗因當因時而異，因專業而異，緊跟時展，做到充分吸引學生的注意力，激發學生的興趣，調動學生積極性，讓學生積極主動的去想辦法解決問題，而非被動的接受。這樣的數值實驗才能有效培養學生應用能力。
數值分析的內容廣泛，包含插值擬合、數值積分與數值微分、數值代數、微分方程數值解法、非線性方程與方程組的數值解法，最優化等等。下面以插值擬合、微分方程數值解法、非線性方程數值解法、最優化的相關實際問題為例，研究數值實驗的設計。
2 結束語
通過以上案例發現：如果知道解決實際問題的數值分析方法，利用MATLAB能很方便的求出實際問題的結果。但是對于多數學生而言，由實際問題無法得到數學模型，更無法知道相應的數值分析方法。因此在學生有一定數值分析基礎后，才能引入這些實際問題的數值實驗；同時選取的數值實驗必須為學生深入研究預留了空間，因為在講解相關的知識背景，詳細分析問題，建立模型后，對學生進行的是分層次指導解決這些問題（即或利用MATLAB函數命令簡單編程計算，或利用MATLAB的專用工具箱計算，或設計算法流程利用MATLAB編程計算）；這樣每一個學生都能參與到實驗中來，每一個學生在實驗中都有收獲；最后撰寫實驗報告，闡明實驗的目的、要求、過程、收獲等.如此，通過這些數值實驗，學生既自己動手解決實際問題，培養了其綜合應用能力，又讓其充分體會到數值分析的魅力，為大學生數學建模競賽和后續專業學習打下堅實基礎.當然，隨著技術的進步，不同專業的數值分析應用實驗需要教師不停的探索和研究。
參考文獻：
[1] 徐翠薇，孫繩武.計算方法引論[M]. 3版.北京：高等教育出版社，2007.
[2] 呂國英.算法設計與分析[M]. 2版.北京：清華大學出版社，2009.
[3] 張德豐等. MATLAB數值分析[M]. 2版.北京：機械工業出版社，2012.
[4] 杜廷松.關于《數值分析》課程教學改革研究的綜述和思考[J]. 大學數學，2007，23（2）：8-15.
[5] 張濤，. 數值實驗在《數值分析》教學中的重要性分析[J]. 長江大學學報：自然科學版，2009，6（1）：371-372.
1計算方法
1．1槳葉／機身非定常面元
1．1．1面元基本原理除物面附近及尾流區外，旋翼流場可假設為無粘、無旋、不可壓。在慣性坐標系下，連續方程可表示成速度勢的函數［17］，即式中SB與SW分別為物面（槳葉或機身）和尾跡渦面，n為物面外法線單位矢量，r＝（x，y，z）為空間點位置。
1．1．2邊界條件物面邊界條件要求相對于物面的法向速度為0，遠場邊界條件要求物體對流體的擾動在無限遠處為0，即假設物體表面由N個面元組成，尾跡渦面由Nw個面元組成，采用等強度四邊形偶極子面元，則式（2）可表示成如下：
1．1．3面元壓力旋翼流場確定之后，可根據非定常Bernoulli方程，通過速度勢和物面速度計算壓力分布。非定常項／t可通過求解物體表面速度勢得到。對于機身，非定項主要來源于槳葉和旋翼尾跡的影響。槳葉影響可通過速度勢直接求解，而尾跡影響為尾跡對面元的誘導速度與尾跡自身速度之積［19］。
1．2時間步進自由尾跡為求解槳葉和機身面元強度分布，在解式（5）或式（8）之前需計算旋翼尾跡。本文采用時間步進自由尾跡［12－14］。時間步進自由尾跡基于不可壓假設，并把旋翼尾跡漩渦簡化為直線渦線。旋翼渦量場可由三維不可壓粘性Navier－Stokes方程描述，表示成速度－渦量采用有限差分近似時間和空間導數求解式（14）。渦線位置由時間步進格式求解得到。文中采用二階精度的預估－修正格式（PC2B）［12－14］。
1．3旋翼槳葉運動方程旋翼尾跡和槳葉面元匯／偶極子分布與槳葉的揮舞運動方程緊密相連，因此在描述旋翼尾跡時需求解槳葉的揮舞運動。根據槳葉揮舞鉸力矩為零建立剛性槳葉揮舞運動方程。槳葉揮舞運動可表示成一組常微分方程，并采用四階Runge－Kutta求解［14］。
1．4槳葉面元／尾跡耦合為計算旋翼尾跡的畸變效應，采用全展渦線代替偶極子面元尾跡。旋翼槳葉由非定常面元構成，槳葉脫出的尾隨渦由尾隨偶極子面元構成，旋翼尾跡則由連接于尾隨渦的全展渦線構成，并從槳葉尾隨偶極子面元中脫出（如圖1）。基于槳葉后緣Kutta條件及尾跡偶極子面元強度與渦線渦量強度等價原則，建立槳葉面元與尾跡之間的聯系。偶極子面元與渦線等價原則可表示成如下。在各時間步，旋翼尾跡渦線強度由槳葉面元強度決定，同時槳葉面元的匯／偶極子強度又與旋翼尾跡渦線有關，由此確保槳葉非定常面元與旋翼尾跡的緊密耦合。
1．5旋翼尾跡／機身干擾低速前飛狀態下，機身浸潤在旋翼尾跡中，因此旋翼尾跡渦線將靠近機身表面。由于機身的阻塞效應，旋翼尾跡渦線靠近機身表面的速度減小，而切線速度增加，此時機身非定常壓力主要來源于旋翼尾跡。由式（11）可知，非定常項由尾跡渦線移動速度和尾跡誘導速度構成，因此旋翼尾跡幾何特性對旋翼尾跡／機身干擾影響顯著。由于機身表面載荷與旋翼尾跡幾何密切相關，因此旋翼尾跡渦線靠近機身表面的運動特性就顯得非常重要［22］。為滿足機身表面無穿透條件，并模擬渦線靠近機身表面的加速現象，文中采用渦線鏡面法。與二維點渦鏡面類似［22］，尾跡渦線由兩點直線構成，因此可通過渦線中點的矢量鏡面得到鏡像渦線，鏡像渦線渦量為Γ′＝－Γ（如圖2）。在各時間步，通過槳葉和機身非定常面元同步求解，得到槳葉和機身的非定常氣動力，而后推進旋翼尾跡，由此計算旋翼／機身非定常氣動干擾。
2計算結果與分析
為驗證本文旋翼／機身非定常氣動干擾分析方法的準確性，文中將計算前飛狀態的Maryland、ROB－IN旋翼／機身干擾下的機身非定常壓力分布，并與可得到的實驗值、CFD計算結果對比驗證。隨后分析前飛速度、旋翼與機身高度對非定常氣動干擾的影響。
2．1Maryland旋翼／機身干擾本算例為前飛狀態下的Maryland旋翼／機身干擾試驗［23］，旋翼系統由4片直徑為1．65m的矩形槳葉鉸接構成，槳葉線性負扭為－12°，翼型為NASARC310和RC410，弦長為0．0635m，旋翼轉速為1860rpm，機身長度為1．94m，機身最大截面直徑為0．254m，機身尾梁與機身截面直徑之比為1∶2．5。槳轂中心與機身重心高度為0．24m。機身壓力傳感器分布如圖。旋翼／機身干擾下的各傳感器非定壓力隨槳葉方位角變化歷程如圖4。從圖4中可以看出，本文計算方法計算得到Maryland旋翼／機身干擾下的非定常壓力時間變化歷程與實驗測量值吻合較好。圖4中各傳感器非定常壓力隨方位角的變化均表現為4Ω周期波動，此倍頻與旋翼槳葉片數相同，由此說明槳葉通過機身上方所產生的顯著非定常干擾效應。機身頭部（傳感器1）非定常壓力呈現出類正弦波動，主要原因為傳感器1在旋翼下方，受槳葉通過性影響顯著。傳感器9、10在旋翼下方之外的尾梁，受到槳葉通過性影響減小，而主要受到旋翼尾跡與尾梁干擾影響，因此非定常壓力呈現鋸齒形狀。傳感器9比傳感器10更靠近旋翼，因此傳感器9的非定常壓力受到槳葉通過性影響更顯著，表現的類正弦特性更顯著。傳感器11、12在尾梁左右兩側，主要受到旋翼尾跡／機身干擾影響，表現出鋸齒形狀。從圖4（b）、（c）中可以看出，傳感器9的壓力峰值相位超前于傳感器10，主要原因為尾梁傳感器9比傳感器10更靠前，旋翼尾跡將先靠近傳感器9。但隨著誘導速度的向下作用，尾跡距傳感器10的距離更小，因此負壓峰值更大，旋翼尾跡／尾梁干擾更顯著。從圖4（d）、（e）中可以看出，尾梁左側傳感12的非定常負壓峰值大于尾梁右側傳感器11，主要原因為旋翼右旋，旋翼尾跡貼近尾梁的左側，因此對傳感器12的干擾作用大于右側的傳感器11。
2．2ROBIN旋翼／機身干擾本算例為前飛狀態下的ROBIN旋翼機身干擾試驗［24］。旋翼系統為2MRTS（2－MeterRotorTestSystem）縮比旋翼［25］，機身為流線型機身。2MRTS旋翼由4片矩形槳葉鉸接構成，槳葉半徑為0．861m，弦長為0．0663m，線性負扭為－8．0°，翼型為NACA0012翼型，旋翼轉速為2000r／min，前進比為μ＝0．151。機身長度為1．999m，槳轂與機身重心垂直距離為0．322m，旋翼周期變距由風洞試驗據得到［25］。各片槳葉由弦向60段和展向20段面元組成，旋翼系統共由4800個面元構成，ROBIN機身由10842個面元組成。機身頭部、發動機艙、尾梁、機身左右兩側壓力傳感器分布如圖5所示。ROBIN旋翼／機身干擾下的旋翼尾跡如圖6所示，從圖中可以看出，旋翼左右兩邊形成比較明顯的槳尖渦。由于機身的排斥作用，旋翼中間尾跡向上、左右兩側移動，但尾梁后段仍然浸入在旋翼尾跡中，因此將產生顯著的旋翼／機身干擾。槳葉／機身壓力分布如圖7所示，機身頭部和發動機艙前、后緣部分產生較大壓力。前飛狀態下，由于旋翼前行邊和后行邊槳葉相對來流的非對稱，因此需通過周期變距改變槳葉的槳距以保證整機左右平衡，并由此導致前行邊和后行邊槳葉氣動環境、槳葉脫出渦量不一致，從而引起旋翼尾跡的非對稱。由于旋翼尾跡非對稱和槳葉位置的變化，導致機身前后、左右兩側壓強非對稱，由此產生時變載荷。ROBIN機頭頂部、發動機頂部、尾梁頂部、機身兩側非定常壓力隨槳葉方位角變化如圖8。從圖8中可以看出，本文計算方法計算得到ROBIN旋翼／機身干擾下的各傳感器非定常壓力間變化歷程與實驗值［24］和CFD計算結果［4，26］比較吻合，由此驗證本文計算方法的可靠性。由于旋翼系統采用4片槳葉，因此圖8中各傳感器非定常壓力隨方位角變化歷程均呈現4Ω的周期特性。從圖5可以看出，各傳感器均位于旋翼下方，因此各傳感器非定常壓力主要表現為槳葉通過性影響。為反映機身各處壓力的變化特性，機身各傳感器非定常壓力峰值相位和幅值如表1．從表1可以看出，機身頭部傳感器6非定常壓力峰值出現在槳葉通過機身后，而發動機艙傳感器22和尾梁頂部傳感器15非定常壓力峰值出現在槳葉未通過機身前，主要原因為機身頭部距離槳尖平面距離更大，阻塞效應較小，且存在阻塞滯后，并由此導致機身頭部非定壓力幅值小于尾梁頂部。由于傳感器22處于槳根下方（圖5），因此受槳葉通過性的影響小于傳感器15。由于槳葉右旋轉，導致機身左側流場阻塞，機身右側流場擴展，由此導致左側傳感器13的壓力幅值大于右側傳感器19，且右側峰值相位滯后與左側。由于機身頭部、尾梁、機身左、右側非定常壓力幅值與相位的差異導致機身力與力矩的非對稱，由此產生4Ω周期激勵載荷。
2．3前飛速度對旋翼／機身干擾影響以ROBIN旋翼／機身干擾為基本算例，分析前飛速度對機身非定常壓力的影響。從圖9中可以看出，隨著前飛速度的增加，旋翼載荷增加，槳葉通過性對機身非定常壓力影響增加，由此導致機身頭部、發動機艙、機身右側壓力幅值均增加，且隨著前飛速度增加，壓力幅值增加速率增加。但尾梁傳感器15壓力幅值隨前飛速度的增加而先加后減小，原因為前飛速度較小時，尾跡／尾梁干擾顯著，速度增加導致尾跡對尾梁的誘導非定常項影響增加，因此壓力幅值增加，但隨前飛速度的繼續增加，尾跡距尾梁的距離增加，因此對尾梁的影響減小。機身與發動艙連接處傳感器19的壓力幅值隨前飛速度增加而減小，原因為傳感器在前行槳葉下方，前飛速度增加，為滿足配平條件，需減小前行槳葉槳距，前行槳葉槳根載荷減小，因此旋翼槳葉通過性影響減小。
2．4旋翼與機身距離對旋翼／機身干擾影響以ROBIN旋翼在前進比為0．15狀態下為基本算例，旋翼與機身高度分別增加5％、10％、15％、20％、30％后機身各部分壓力幅值的變化如下。從圖11中可以看出，隨著旋翼與機身距離的增加，槳葉通過性對機身非定常壓力影響減弱，由此導致機身頭部、尾梁頂部、機身左右側各處壓力幅值均減小。隨著旋翼與機身距離的增加，壓力幅值減小速率逐漸減小，距離增加20％，傳感器幅值減小為參考值的80％以下。但發動機艙頂部傳感器幅值隨旋翼與機身距離的增加而先減小后增加，主要原因為旋翼與機身距離的增加，槳葉通過性影響減小，因此壓力幅值先減小；旋翼與機身距離的繼續增加，方位角為270°處槳葉的槳尖渦將貼近發動機艙頂部，由此導致壓力幅值增加。
3結論
關鍵詞: 數值分析 數學建模 Matlab
數值分析又稱計算方法,是一門與計算機使用密切結合的實用性很強的一門課程,重點研究如何運用數值計算方法去處理實際工程問題,因此數值分析在科學研究、工程建設和經濟建設等很多方面有著廣泛的應用。在信息科學和計算機技術飛速發展的今天,這門課程中的數值方法更顯得極其重要,但是對多數學校來說,還沒有引起對這門課足夠的重視,而且在數值分析的教學過程中都存在很多不足。不少學者也討論過我國高校中數值分析課程的教學情況,其中存在一些普遍問題,例如學生理論學習模式化、實踐能力不夠、缺乏應用性,學習過程中學生感覺到枯燥或者學習效果不佳,學校軟、硬件設施無法滿足學生的上機實習等。如何更好地開展這門課程的教學工作,對于我們來說是一個巨大的挑戰。下面我們來談談在教學過程中遇到的幾個問題。
1.理論基礎知識扎實,同時采用啟發式教學
課程中的很多公式是推導出來的,推導過程比較煩瑣,得到的公式也比較冗長,而且比較難記,對于已經復雜并且很冗長的數值公式,還需要進一步進行抽象的理論分析,包括算法的收斂性如何,數值算法是否穩定并進行誤差分析,以及分析算法的空間和時間復雜性等,同時還涉及如微積分、線性代數、常微分方程等。過多地強調數學理論證明,大多數的學生覺得這門課很難,學得很枯燥,也感覺不到樂趣,從而越來越厭煩學習這門課程。
因此,我們要將“因材施教”的理念落到實處。方法的講授應該盡量地從實例中提出問題,引導學生去思考如何運用數學知識去構造解決的方法,然后給出相應的數學理論。并且,給出一種方法,可以換位思考,激發學生思考是否能用另外的已學方法來求解。這樣不僅能復習已學的知識,而且能鞏固各種知識之間的聯系,還可以啟發學生把學過的知識學以致用,真正了解學習帶來的樂趣。
2.將數學建模的思想融入到教學過程中
數值分析是對實際問題的數值模擬方法的設計、分析與軟件實現的理論基礎。要解決具體的實際問題,首先需要建立起適當的數學模型,將實際問題的解決歸結為相應的數學問題的求解,然后對所歸結的數學問題建立相應的數值方法。這樣就可以以實例啟發學生弄清為什么要進行數值分析、應該如何引進數值方法進行分析,建立一種數值分析的方法后,哪些問題是值得且必須研究的。例如在汽車、飛機等的外形設計過程中,利用樣條技術設計的外形越來越光滑、美觀。學生了解了樣條插值的實際應用背景后就會對樣條插值的理論更感興趣,也會更有動力來學。
將數學建模的思想融入到數值分析教學過程中,要求我們必須有一個合適的切入點,不能用數學建模課的內容過多占有數值分析課的教學,因此精選只涉及相應數值分析理論和方法而又能體現數學建模思想的內容,既能吸引學生又是學生以后可能碰到的案例,將其融入到數值分析課程中是十分重要的。下面具體舉兩個例子,插值方法可以引入人口增長的模型和設計公路平面曲線的問題,常微分方程的差分方法可以引入導彈追蹤和估計水塔的流量問題,方程求根的迭代法可以引入一般戰爭模型,線性方程組的解法可以引入投入產出模型和小行星軌道問題等。
3.結合Matlab進行實踐教學
在結合多媒體教學的過程中,盡量地在講解數學模型的過程中,無論是問題的引入還是算法的講解和實現,以及結果盡可能地轉化成圖形等一些可視的結果展示給學生,以激發學生的學習興趣,引人入勝,Matlab軟件的可視化功能能夠實現這一點。
在計算機技術飛速發達的今天,只要有效地把教學過程和相關的計算機技術結合起來,就能夠做到減輕教師教和學生學的負擔,優化學習環境,實現高效教學。在一些數值分析教材中一些常用的算法都已經有了現成的程序,因此在授課的過程中,對這些算法進行展示時,要讓學生從中學會如何將一個算法轉變成一段程序。鼓勵學生自己根據算法寫出程序流程圖,然后使用Matlab語言將其轉變成程序,將自己所得程序與課本中的結果進行比較分析,這個過程有助于學生更好地理解算法,增強學生動手實踐的自信心。
4.結語
數值分析是研究數學模型的數值計算方法。隨著電子計算機的迅速發展、普及,以及新型數值軟件的不斷開發,數值分析的理論和方法無論是在高科技領域還是在傳統學科領域,其作用和影響都越來越大,實際上它已成為科學工作者和工程技術人員必備的知識和工具。
對于理工科的本科學生而言,它的理論和實踐知識對學生的要求都比較高。因此要讓學生學好這門課程,需要在教學中采用一些技巧性的教學方法,比如采用啟發式的教學方法,融入數學建模的思想,以及結合Matlab進行實踐教學等。這樣可以調動學生主動學習的積極性,提高學生的綜合素質,使學生真正學好這門課程。
參考文獻:
[1]趙景軍,吳勃英.關于數值分析教學的幾點探討[J].大學數學,2005,21(3):28-30.
[2]孫亮.數值分析方法課程的特點與思想[J].工科數學,2002,18(1):84-86.
Abstract： Based on the bubble dynamic equation under the consideration of liquid surface tension, viscosity and radiative resistance, this essay adopted numerical simulations to investigate single cavitation bubble dynamics with different kinds of acoustic driving frequencies and bubble initial radiuses.
關鍵詞： 超聲空化；空化氣泡；數值分析
Key words： ultrasonic cavitation；cavitation bubble；numerical analysis
中圖分類號：G30文獻標識碼：A文章編號：1006-4311（2011）01-0196-02
0引言
隨著科學技術的發展，超聲已在眾多領域得到了廣泛的應用，在這些應用中，超聲空化是引發各種物理、化學和生物效應的主要機理，這些效應與瞬態空化氣泡崩潰時所產生的高溫高壓等現象有關。在研究單一空化氣泡動力學過程的方法中，數值分析是除理論和實驗方法之外的一種研究方法，至少有兩方面原因表明它是必要的。首先，由于氣泡運動過程中高度的非線性，使得從理論上建立能夠精確描述空化過程的方程實際上是不可能的，其次，微米級大小的空化氣泡半徑和持續時間為微秒至納秒級的氣泡運動周期使得實驗測量也難以進行。本文基于考慮了液體表面張力、液體粘滯性和輻射阻尼的氣泡運動方程，采用數值分析中Runge-Kutta方法研究在不同聲場頻率和氣泡初始半徑條件下單一空化氣泡的運動過程。
1氣泡動態的數值分析
考慮了液體表面張力、液體粘滯性和輻射阻尼的單一空化氣泡運動方程：
R+=P+-P-P-
-+P+-P（1）
方程（1）是二階常微分方程，使用Runge-Kutta方法求解時應當用替換法化為形如y′=f（x，y）的一階微分方程組，再加以使用。方程（1）可化為下列一階方程組：
y=R′y′=f（t，R，R′）=-（R′）+P+-P-P--R′+P+-PR?佐=R，y?佐=R′?佐=0
若設時間步長為h，則使用Runge-Kutta方法求每個離散時間點上對應氣泡半徑大小Rn的遞推公式為：
R=R+hR+k+k+kR=R+k+2k+2k+kk=hft，R，Rk=hft+，R+R，R+k=hft+，R+R+k，R+k=hft+h，R+hR+k，R+k
設聲場激勵波形為P=Psin（ωt），若聲壓P
f=P+-（2）
同時假設氣泡運動過程為等溫過程，即方程（1）中的泡內氣體多方指數n=1。計算時與液體相關的各參數取值分別為：液體密度ρ=1000kg/m3，液體表面張力系數σ=0.076N/m，液體粘滯系數μ=0.001kg/（m•s），液體中聲速c=1481m/s，液體中的靜壓P=1.013×105Pa。
下面研究單頻聲場激勵下，聲場頻率f和氣泡初始半徑R對氣泡動態的影響。具體為計算f和R在不同的取值條件下，空化氣泡半徑隨時間的變化曲線，即R（t）曲線。
1.1 聲場頻率對氣泡動態的影響給定氣泡初始半徑R，通過改變聲場頻率f的大小，討論fa的變化對氣泡動態的影響。設聲場激勵為P=-Psin（2πft），P=5.0×105Pa，R=0.6μm，由公式（2）得到其自然共振頻率f=7.55MHz。若取f=20MHz，則反映氣泡動態的R（t）曲線為圖1所示，氣泡在多個聲波周期內做復雜振蕩。
若取f=7MHz，則R（t）曲線為圖2所示，氣泡在一個聲周期內即趨向崩潰。需要指出，這里所說的趨向崩潰是指氣泡半徑在增大到最大值后急劇向R=0趨近。
氣泡在不同聲場頻率激勵下其趨向崩潰的程度是不同的。圖3顯示了聲場頻率從4MHz變化到10MHz時，氣泡在趨向崩潰時的半徑與其初始半徑之比R/R0的變化趨勢。
通常認為:當超聲波頻率與氣泡的自然共振頻率相等時，超聲波與氣泡之間就達到了最有效的能量耦合，氣泡將迅速崩潰。但數值計算的結果表明，當f=f時，氣泡半徑在通常所指的崩潰階段趨向0，但不為0；當f小于f至一定限度時，氣泡半徑才在10-5數量級上趨向0，這時可以認為氣泡徹底崩潰；當f大于f至一定限度時，氣泡可在多個聲波周期內穩定振蕩，且振蕩波形復雜無規律。
同時大量實驗研究表明，隨著頻率升高，聲空化過程變得難以發生。對這種現象的定性解釋為：頻率增高，聲波膨脹相的時間相應變短（如f=20kHz，其膨脹時間為25μs；如f=20MHz，其膨脹時間為25ns），氣泡核來不及增長到可產生效應的空化氣泡，或者即便空化氣泡可以形成，但由于壓縮相時間也短，空化氣泡可能來不及收縮至發生崩潰。為使在較高頻率下產生空化，可以提高聲強，即空化閾值將隨頻率升高而增大。此外，從聲波的傳播特性可知，頻率升高，聲波的傳播衰減將增大，這也使得空化強度減弱以及可能發生空化的區域減小。
1.2 氣泡初始半徑對氣泡動態的影響給定聲場頻率f，計算氣泡取不同初始半徑R時的動態曲線R（t）。設P=P=1.013×105Pa，f=20KHz，這是超聲工業清洗及聲化學中較常使用的頻率，按照公式（2）與其對應的共振氣泡半徑為R=150μm。圖4為R從60μm變化到160μm氣泡趨向崩潰時R/R的變化趨勢。同樣可以看出，當聲場頻率一定時，在該頻率上自然共振的空化氣泡并沒有徹底崩潰，而是那些半徑小于自然共振半徑至一定限度的氣泡才趨向徹底崩潰，半徑大于該自然共振半徑的氣泡將持續振蕩若干周期而不崩潰。
一般情況下，當時，計算得出下述結論：對初始半徑大于共振半徑的氣泡，將發生復雜的持續振蕩，一般不會崩潰；對初始半徑小于共振半徑的氣泡，隨著聲壓負壓相的到來而不斷增大，當聲壓正壓相到來時，氣泡先因慣性繼續生長到最大半徑，然后迅速收縮，直到崩潰。
2結論
本文采用Runge-Kutta數值分析方法研究在不同聲場頻率和氣泡初始半徑條件下單一空化氣泡的運動過程，數值分析結果表明，當給定氣泡初始半徑大小時，聲場頻率在小于氣泡自然共振頻率至一定限度時，氣泡將迅速崩潰，而大于該共振頻率時，氣泡將持續振蕩而不崩潰，即隨著聲場頻率升高，聲空化將難以發生；當給定聲場頻率時，只有其半徑小于與該頻率對應的氣泡自然共振半徑至一定限度的氣泡才會徹底崩潰，半徑大于該自然共振半徑的氣泡將做持續振蕩。
參考文獻：
[1]馮若.超聲手冊[M].南京：南京大學出版社，2001.
[2]錢夢J，程茜，葛曹燕.單泡聲致發光中氣泡的動力學特征――振子模型[J].聲學學報，2002，27（4）：289-294．
[3]李信真，車剛明，歐陽潔，封建湖.計算方法[M].西安：西北工業大學出版社，2000.
關鍵詞地基承載力，塑性力學上限，最優化方法。
1前 言
地基承載力、土壓力和邊坡穩定是土力學的3個重要領域。這 3個問題都基于共同的極限平衡分析原理 ,可以采用相同的分析方法。但是 ,在長期的實踐中 ,這 3個領域各形成了自己的體系。在地基承載力領域 ,目前常用的計算方法仍然是基于 Prandtl解的各種經驗修正公式。應用塑性力學上下限原理 ,在建立地基承載力、土壓力和邊坡穩定分析統一的理論和方法方面作了大量的工作 ,但其有關的研究一直是在變分原理基礎上進行的 ,因此 ,難以擴展到具有復雜邊界和分層土體的實際工程問題中。曾提出一個基于滑楔破壞模式的分析方法 ,其普遍適用性還有待進一步論證。顯然 ,只有開發數值分析的方法 ,方可使大部分實際問題方便地獲得解答。
近期 ,筆者在二維領域應用塑性力學上限定理進行邊坡穩定的理論研究[4 ,5]。該方法從變形協調出發 ,對于一個設定的滑裂面和斜分條模式 ,建立協調的速度場 ,根據外力功和內能耗散相平衡的原理確定相應的安全系數或加載系數 ,然后應用最優化方法 ,確定對應于最小安全系數的那個臨界滑裂面和斜分條模式 (以下簡稱能量法 )。這一方法在精確地確定邊坡穩定安全系數方面獲得了成功。由于地基實際上是一個坡度為零的邊坡 ,將該成果推廣到地基承載力 ,是一個十分具有理論和實用價值的課題。
2極限分析法的理論基礎和計算步驟
2.1上限定理的基本命題
在邊坡穩定和地基承載力分析領域 ,對上限定理的描述可以用下面的命題表達(圖 1):
在塑性區 Ω*,給出一個機動可能的應變場εij* ,并在滑裂面 Γ* 上給出一個相應的速度場 V*,那么,按照下式計算獲得的外荷載T* 將比一個包含有真實的塑性區 Ω和真實的滑裂面 Γ的臨界荷載 T大或與其相等。
∫Ωσi*jεi*j dΩ +∫vdDs* = WV* + T* V*(1)
上式左邊的第一、第二項分別為塑性區內和滑裂面的內能耗散；W為塑性區土體重。因此 ,在諸多協調的位移場中給出最小的 T* 的那個一定最接近真實的臨界荷載 T。
在地基承載力問題中 ,通常定義加載系數 η* 為
η* = (T * - T0)/ T0 (2)
其中T0為地基的實際承受的外荷載 ,那么上限圖定理的命題具體化為尋找一個使 η* 獲得最小值 η的應變場和速度場。
2．2計算內能耗散
如果材料遵守莫爾-庫侖破壞準則和相關聯的流動法則，則可確認沿滑面的速度V與滑面夾角為摩擦角φ，單位面積內能耗散為（圖2）：
d D = (c cosφ u sinφ)v(3)
其中c為凝聚力；u為孔隙壓力；V為滑塊沿滑面的在單位荷載增量下產生的相對位移 ,通常稱變形速率。
2.3計算多塊體破壞模式協調的速度場
對某一邊坡的塑性區 ,將其用一系列傾斜的線分成若干楔塊 ,每一楔塊都視為剛體 ,其變形速率為 V。圖 3示出 3個塊體的系統。 V與滑面夾角為 φ,與右邊相鄰塊體的相對速度為 Vj ,V j與該兩塊體交界面的夾角為 φj。內能耗散發生于該楔塊的底面和楔塊間的界面 ,在剛體內為零。
根據位移協調要求 ,可以得到
1 (4)
1 (5)
其中Vl 和 Vr 分別為左側和右側條塊的速度 ;θj =π/2-δ+ φj, θl=π+αl-φl,θr =π+αr -φr;α為底面與 x軸正向夾角 ;δ為側面與 y軸正向夾角 ;θ為速度與正 x軸的夾角。如果將條塊的寬度取為無限小 (圖 4),還可通過積分獲得滑面上坐標為 x的條塊的絕對速度和相對速度。
V0為左端點（x=x0）的速度。在滑裂面上第k個α或φ 發生突變。上標l和代表該不連續點左和右的物理量。計算從第一個界面開始，到第 n-1個界面終止。這樣，對滑面上橫坐標為x的任意一點，其條塊絕對速度V和條塊側面的相對速度Vj 都可表達為滑面左端點x=x0處的速度V0的函數。將獲得各條塊的絕對速度和相對速度代入式 (3)再代入式 (1),其中式 (1)左側第一項可通過將Vj替代式（3）中 V獲得。消去左右側 V0 ,就可求解加載系數 η*如下。
定義
其中d W =土條重量 ;T0x, T0y分別為 T0在 x和 y軸的分量 ;L為土條側面長度 ;η′為水平地震力系數。式 (10)最后一項計及滑面上 (n -1)個 α或φ的不連續點相應的界面上的內能耗散。由式 (2)定義的加載系數可通過下式計算 :
η* = G/Gb(12)
2.4求解臨界滑動模式
陳祖煜和邵長明曾詳細介紹了對傳統的極限平衡法計算最小安全系數和臨界滑裂面的數值分析方法。最優化方法為使用計算機搜索臨界滑動模式創造了條件 ,這些研究成果可以方便地推廣到本文介紹的極限分析方法中。所不同的是 ,滑動模式和垂直條分法相比 ,增加了一個土條界面傾角 δ。每一條塊的 δ也將成為自由度。最優化方法將最終找到相應最小加載系數的滑裂面和斜分條模式。具體計算步驟通過下節 [例 2 ]介紹。
3驗證
為了驗證上述推導的正確性 ,下面通過兩個例題進行分析探討。
[例 1 ]對具有垂直表面荷載的例題 (圖 5),索科洛夫斯基 ( Sokolovski , 1954)給出的臨界垂直荷載 q的計算公式如下 :
的 q=c ctgφexp-2x)tgφ-1
其中 χ為邊坡斜面相對水平線的夾角。相應的臨界滑裂面由三段組成 ,AB , CD為直線 ,分別與邊坡線和坡頂線夾角為 μ。
μ=
BC為一對數螺旋線 ,其左右邊界線 BO和 CO分別與邊坡線和 y軸線夾角為 μ。
當邊坡處于極限狀態時 ,加載系數 η=0。對 AB ,BC和 CD段分別進行積分 ,按式 (13)確
定的 q將使按式 (10)確定的 G為零。
這一實例說明 ,本文提出的上限定理的命題可通過解析解獲得印證。
[例 2 ]某一坡角為 35°的均質邊坡 ,其水平頂面上作用一均布荷載 ,荷載方向相對鉛直
線夾角為 δ′(圖 6)。根據索科洛夫斯基 (1954)提出的滑移線方法 ,此題的理論破壞面由直線 AB , CD和對數螺旋線 BC組成 , CD和 CO分別相對鉛直線夾角 μ+ρ和 -(μ+ρ)。其中 μ = 45°-φ/ 2 ;ρ為大主應力相對鉛直線的夾角。
主要參數 :c = 720kPa ,φ = 37°,χ = 35°,δ′= 24°,理論解提供的解答是 q =6 . 228 MPa , ρ= 28. 4°。理論的滑裂面和土條側面示于圖 6線 4。滑裂面通過聯結 4個點的樣條函數形成。對設定的初始滑裂面 1和相應的斜分條模式使用式 (12)求得 η3 =0127。從這個滑裂面開始 ,進行最優化方法計算最終得臨界滑裂面和條間界面 (滑裂面 3 ,虛線 ),相應 η=01019。滑裂面 2是優化計算過程中通過隨機搜索獲得的滑裂面。如果用 5個點來模擬該滑裂面 ,則可得到 η=010028。與理論解相比 ,無論是最小加載系數 ,還是臨界滑裂面和臨界條間界面均十分接近。
通過 [例 2 ]說明 ,應用最優化方法可以自動找到相應最小加載系數 η的臨界滑裂面和相應的斜分條模式。
4能量法在地基極限承載力計算中的推廣
4．1傳統的承載力計算方法
地基極限承載力的計算包括兩個方面，一方面是允許位移的校核 ,另一方面是極限承載力的計算。對于后者， Prandtl于 1920年根據塑性力學理論導出了剛性基礎壓入無重量土中滑裂面的形狀及其相應的極限承載力計算公式。由于數學上的嚴格解答在大部分的實際問題中是不可能得到的 , Terzaghi ,Meyerhof ,Vesic等眾多學者在 Prandtl解的基礎上對承載力理論進行了研究和發展 ,最終形成地基極限承載力的近似解答。這一近似解答的一般表達式為
Nγ為一半經驗數據 ,可從地基規范承載力表中查取或用半經驗公式 (表 1)計算 ;B為基礎寬度 ;D為基礎埋深 ;γ為土容量 ; qu為地基的極限承載力 ,即 T在單位寬度上的強度。
4.2能量法在地基極限承載力計算中的應用
選取寬度 B = 17m的條形基礎進行計算分析 ,相應參數為 :c = 144. 5kPa ,γ=0. 0kN/ m3。這個例子針對土的不同內摩擦角 φ值進行計算。對于具有理論解的實例 ,使用式 (13)獲得的 qu應保證使用式 (12)獲得的 η的最小值為零。圖 7示出 φ =0°和φ = 20°兩種情況。使用同樣的初始破壞模式如圖 7 (a),應用最優化方法獲得的臨界破壞模式分別如圖 7 ( b)和 7 (c)所示 ,η分別為 01004和 01008。計算機在搜索最小 η值時 ,準確地將中部的土條側面收斂于地基的左側點 ,由此將滑裂面分為 3個區域 :荷載作用下面的三角形區域 ,對應于理論上的主動 Rankine區 ;條間界面一端收斂于一點所對應的放射形區域 ,對應于理論上的 Prandtl區;放射形區域另一端的另一個三角形區域對應理論上的被動 Rankine區。 φ=0°時 ,滑裂面形狀接近于圓弧 ;φ≠0°時 ,滑裂面形狀為兩段直線接一段對數螺旋線。這就說明 ,采用建立在塑性力學上限 解基礎上的地基極限承載力數值分析方法直接獲得了理論上嚴格的地基極限承載力解答。
表 2將一系列φ值的計算成果與理論解對比，可見成果的準確性相當穩定。所得的η值與理論值的誤差均在1%左右，而且自動搜索得到的臨界滑裂面形狀也與理論解一致。
表2 上限數值解和理論解成果對比
注：qn和 qp分別為根據數值解和理論解獲得的極限承載力 ;η= ( qn/ qp -1);q的單位為 kN/ m。
表 3比較了容重不為零情況下表 1所列的各種經驗方法的準確性 ,并示于圖 8。可見 W1F Chen的公式計算結果與采用 Prandtl的經驗公式求得的結果最為接近 ,但是在 φ值超過 25°后 ,η值為負 ,意味著高估了地基的承載力。而 Mayerhof以及 Terzaghi的方法則偏于保守。
表 3對 γ≠0情況各種不同的經驗公式和數值解成果對比
注 :下標 n,v, m,t,c分別代表數值解、采用VesicMeyerhofTerzaghi 和 W. F. Chen方法的計算成果，q的單位為 kN/ m。，
在有容重且有埋深的條件下 ,式 (15)的經驗公式將基礎兩側埋置深度以內的土重以連續均勻分布的荷載考慮 ,未能計及這部分土體的抗剪作用 ,因此不可避免地存在著一定的誤差。
[關鍵詞]壓縮式封隔器；擴張式封隔器；封隔器膠筒；有限元分析
1.封隔器的應用現狀
壓縮式封隔器的工作原理：當膠筒承受軸向載荷時，封隔器膠筒將產生徑向大變形，使膠筒與套管之間產生接觸壓力，藉此封隔環空，隔絕產層，以控制產（注）液，保護套管。膠筒所用材料為橡膠，橡膠最顯著的特性就是它的超彈性，即在很小的力的作用下就能產生很大的變形，正是這一性質使得橡膠材料被制成密封元件并在工程中得到廣泛應用。
但是，國內還沒有專門用于深水井下作業的相關封隔工具，因此，研究封隔器膠筒的材料以及研究用于惡劣環境下的井下措施工具以及技術，是國內各大研究單位以及相關高校的當務之急。
2.膠筒材料本構關系
2.1 基本理論
橡膠材料在較短時間內及恒定的環境溫度下通常被視為各向同性不可壓縮材料，其應變能密度函數W是變形張量不變量的函數，即：
（1）
其中：
（2）
（3）
（4）
式中：λ1、λ2和λ3是3個主伸長率。
Rivlin用唯象法從純數學角度出發導出應變能密度函數的最一般的形式為：
（5）
式中：Cijk為Rivlin材料常數。
橡膠材料可以近似認為是不可壓縮的，因此I3=1，式（5）可簡化為：
（6）
取式（6）的前兩項得到Mooney-Rivlin模型：
（7）
其中C10和C01為Rivlin材料常數，均為正定常數。
2.2 Mooney-Rivlin模型材料常數的確定
對兩參數Mooney-Rivlin模型，Mooney-Rivlin材料模型的應力-應變方程可表示為：
（8）
這個方程可用σ/2（λ-λ-2）對1/λ作圖。在1/λ=1時，相應值為C10+C01并且直線的斜率為C01，初始剪切模量與材料常數的關系是G=2（C10+C01）。如果材料被假定為不可壓縮的，那么初始拉伸模量是E=6（C10+C01），G=2（C10+C01）和楊氏模量的值近似等于3倍的剪切模量。
利用橡膠材料的硬度與C10和C01存在一定的關系，根據橡膠材料硬度與彈性模量的實驗數據，經擬合可得：
（9）
利用硬度計可測得橡膠試樣的硬度HA，通過式（9）可算出橡膠彈性模量E，根據公式E=6（C10+C01）可以得到C10+C01的值，通過公式d=（1-2μ）/（C10+C01）可以測出模型參數d的值。利用經驗公式C01=0.5C10求解下面的方程就可得到C10。
E=6（C10+C01），因此E=6（C10+0.5C10），則。
由實驗測得在常溫25℃和高溫150℃下膠筒的邵氏A硬度分別為94和91，通過公式（8）、（9）即可計算出的膠筒的Mooney-Rivlin模型材料參數。
3.封隔器膠筒有限元數值模型的建立
封隔器采用兩膠筒形式，上、下膠筒、中心管和套管結構。上下膠筒力學參數和幾何尺寸完全相同，兩膠筒由隔環隔開。由于封隔器中心管、膠筒、套管以及膠筒所受載荷均為軸對稱分布，故取過軸線的剖面建立有限元計算模型， 膠筒采用四節點四邊形超彈單元PLANE182，中心管、套管和隔環采用四節點四邊形等參數單元PLANE42，用面面接觸單元CONTA171和TARGE169模擬膠筒與套管、中心管和剛性隔環間的接觸，并考慮接觸時的摩擦作用，通過實驗測定膠筒與鋼材之間的摩擦因數為0.3，鋼與鋼之間的摩擦因數為0.1，中心管和套管上下兩端采用固定約束，橡膠筒在上端隔環的作用下壓縮膨脹，與套管接觸，隨著坐封壓力不斷增加，膠筒與套管的接觸壓力也隨之增大，從而起到密封套管內上下環空壓差的作用。
建立膠筒有限元模型，對膠筒進行非線性接觸分析。膠筒具有不可壓縮性，因此取膠筒泊松比為0.4996。
4.數值模擬結果分析
4.1 坐封過程中接觸壓力分布規律
為模擬封隔器坐封過程中膠筒與套管內壁接觸情況，穩定溫度150℃不變，分步向膠筒上端隔環施加軸向荷載，觀察膠筒的變形情況。
由于膠筒與套管內壁接觸壓力和摩擦阻力的存在，2個膠筒所受到的軸向壓縮載荷由上至下依次降低，所以盡管上下膠筒的材料和幾何形狀完全相同，但是下膠筒的變形、與套管的接觸壓力均小于上膠筒，坐封完成以后，上膠筒的接觸壓力比較大，下膠筒的接觸壓力比較小。顯然，隨著封隔器坐封載荷的增加，上膠筒與套管的接觸壓力明顯高于下膠筒的，所以上膠筒在密封套管環空上下壓差中起到關鍵作用。
4.2 坐封過程中軸向載荷與膠筒變形的分布規律
根據ANSYS分析輸出結果可以得到在150℃時不同坐封力下膠筒的軸向壓縮位移和軸向載荷的關系。穩定溫度在150℃時，依次施加10t、12t、13t、15t、17t、18t坐封力時膠筒的軸向壓縮位移分別為23.095mm、24.319mm、24.564mm、24.905mm、25.086mm、25.186mm。
實驗發現，隨著軸向載荷的逐步增大，膠筒軸向壓縮位移也漸漸增大，當軸向載荷小于12t時，膠筒軸向壓縮位移隨軸向載荷增大的梯度比較大，當軸向載荷大于12t時，膠筒軸向壓縮位移增大梯度慢慢減小。因此，當軸向載荷過大時，膠筒軸向壓縮位移的變化不是很大。
5.結論
從數值模擬可以看出，封隔器的兩個膠筒與套管的接觸壓力分布為：下膠筒接觸壓力較小，上膠筒的接觸壓力較大，隨著軸向載荷的增大，軸向壓縮量也增大，同時也可以發現，隨著坐封力的增大，膠筒與套管接觸長度逐漸增加，膠筒外表柱面部份徑向變形受限制，膠筒內表面變形如外表一樣向外鼓，當載荷增加時膠筒被壓扁并在最后壓實。并且隨著軸向載荷的逐步增大，膠筒軸向壓縮位移一開始是漸漸增大，增大的梯度也比較大。當軸向載荷繼續增加時，膠筒軸向壓縮位移增大梯度會慢慢減小。因此，當軸向載荷過大時，膠筒軸向壓縮位移的變化不是很大，為了達到良好的密封，應該計算合適的坐封載荷。
參考文獻
[1]曾憲平.對壓縮式封隔器膠筒變形分析[J].石油鉆采工藝，1983，13（4）：65-68.
[2]楊秀娟，楊恒林.液壓封隔器膠筒座封過程數值分析[J].石油大學學報，2003，05-0084-04.
[3]劉永輝.封隔器膠筒性能有限元分析及結構改進[D]西南石油大學，2006
關鍵詞：混凝土 ;隨機骨料模型; 映射網格;細觀數值模擬
中圖分類號：TU528.01
文獻標識碼：B
文章編號：1008-0422（2008）08-0163-03
1前言
混凝土作為重要的建筑材料已有百余年歷史，人們從不同的結構層次來研究混凝土材料的力學性能。（1）宏觀層次。一般認為特征體積相當于3～4 倍的最大骨料體積。當小于特征體積時，當研究對象大于特征體積時，材料假定為均質。（2）細觀層次。從分子尺度到宏觀尺度之間，其結構單元尺度變化范圍在10-4cm至幾cm。在這個層次上，混凝土被認為是一種由骨料、水泥砂漿及二者之間的界面組成的三相復合材料[1]。(3) 微觀層次。材料的結構單元尺度在原子、分子量級，即從小于10-7cm～10-4cm 。材料由晶體結構及分子結構組成。隨著高性能計算機及數值算法的迅速發展，用數值試驗來模擬混凝土材料破壞過程的研究已成為熱點，相對廉價的細觀數值分析有望彌補或替代部分材料試驗[2]
2骨料尺寸
選取300mm×300mm×300mm的立方體試件，并將它們取為二級配混凝土，其中大石(直徑40mm~80mm)、中石(直徑20mm~40mm)進行數值模擬。按富勒理想級配曲線：
（1）
其中D為篩孔直徑；Dmax為骨料最大粒徑；P為能通過篩孔D的骨料的質量百分比。
取Dmax為80mm，求得大石、中石的質量百分比：
由表1可知，粒徑為40mm~80mm的骨料質量百分比為100%-70.7%=29.3%，粒徑為20mm~40mm的骨料質量百分比為70.7%-50%=20.7%。骨料容重取2800Kg/m3，試件的石料用量為1227Kg/m3。適當選擇球體骨料的個數與直徑，使之滿足級配關系并容易隨機填充，本文取大骨料為大石，個數為27,直徑為62.58mm。小骨料為中石，個數為189,直徑為29.16mm。
3三維有限元網格的劃分
隨機骨料模型的網格剖分，長期以來是三維細觀力學分析的難題之一。分網技術簡單地說分為三種：六面體映射或拉伸網分、四面體全自動網分、混合網分。（1）商用軟件中，TRUEGRID、ANSYS ICEM CFD、AVL_FAME等可以完成六面體映射網格劃分。劃分隨機骨料模型時，主要思想是：先生成均勻的有限元網格覆蓋整個分析區域，根據蒙特卡羅方法和骨料從大到小的投放原則，隨機生成一套虛擬的骨料域。根據單元節點是否全部落在、部分落在和不落在骨料域，判別單元是否屬于骨料單元、界面單元和砂漿單元。這種映射網格技術能夠解決復雜形狀骨料的三維模擬，全部是六面體網格，有限元的單元精度高，但界面區域不夠光滑。在計算上，考慮到界面的厚度很小，要采用較小的單元而在全域上采用一樣大小的單元，導致單元個數太多，計算量較大，在微機上計算難度大。（2）多數的CAE軟件都提供四面體網格全自動劃分功能，配合局部加密技術劃分單元，可以建立隨機骨料網格模型，計算量可以接受，但四面體單元的精度較低，空間填充效率不高，且視覺效果稍差。（3）CUBIT、MARC、HARPOON等軟件提供混合分網功能，產生的網格主要為六面體，邊界復雜的區域用五面體、棱柱體或四面體填充。分網效率高，但邊界區域通常有應力集中，需要特別關注的地方卻用的是低精度網格。
本文提出先分塊映射網分再隨機裝配的方案，可以兼顧單元精度、局部加密和計算量的要求。首先，將混凝土試件看成由三種塊體組合而成：（1）包含大骨料、界面和砂漿的塊體，簡稱大骨料塊。如圖1所示。（2）包含小骨料、界面和砂漿的塊體，簡稱小骨料塊。如圖2所示。（3）只包含砂漿的塊體，簡稱砂漿塊。如圖3所示。分別對這三種塊體用映射網分。
以上是球體骨料的情形，類似地，可作出橢球體骨料的映射網格。如圖4所示。
注意到界面的材料剛度與骨料的剛度有較大的差別，為了減少應力集中，將骨料的網格作了適當的加密。
三種塊體都是正方體，使大骨料塊的邊長為66.7mm，小骨料塊和砂漿塊邊長33.3mm，界面單元厚度為1mm，骨料和界面單元全部位于塊體的內部。大骨料塊的邊長為小骨料塊和砂漿塊邊長的兩倍，表面上的單元尺寸相同，顯然這三種塊體可以無縫、無重疊地填充三維空間，而且將重合的節點合并后，各塊之間形成共節點連接的整體有限元網格，如圖5所示。隨機裝配過程利用了MARC軟件的PYTHON語言前處理接口，調用隨機函數以及單元的復制、移動功能。整體有限元共有252424個六面體單元，258877個節點，其中砂漿單元139968個，如圖6所示。骨料單元82922個，如圖7所示。界面單元28512個，如圖8所示。
4材料參數
Drucker -Prager 屈服準則在MSC.Marc有限元軟件中稱其為Linear Mohr-Coulomb準則，屈服條件為[4]：
為應力張量的第一不變量； 為偏應力張量的第二不變量。
MSC.Marc提供了可模擬開裂和應變軟化的損傷模型。各主要參數如表2所示。
5結果分析
按位移控制加載，鉛垂方向分級加拉伸位移，ε表示混凝土試件在加載方向的總應變，圖10~圖14顯示的是等效開裂應變云圖，相當于顯示出了裂紋分布的區域。
圖10顯示當伸長總應變為ε=0.663×10-4時在大骨料的界面上開始出現裂紋，裂紋方向為垂直于加載方向。圖11顯示當伸長總應變為ε=0.773×10-4時在小骨料的界面上也開始出現裂紋，而且大骨料界面上的裂紋繼續增加，裂紋方向為垂直于加載方向。圖12顯示當伸長總應變為ε=0.8673×10-4時，大骨料界面上的裂紋繼續增多，在小骨料的界面上裂紋沒有明顯增多。
圖13顯示當伸長總應變為ε=0.970×10-4時，在小骨料的界面上裂紋有明顯增多。圖14顯示當伸長總應變為ε=1.07×10-4時，在大、小骨料的界面上裂紋足夠多了。雖然裂紋之間還沒有貫穿相連，但材料非線性的求解得不到下一載荷步的收斂解答，可認為材料已達到單向拉伸的極限載荷。
6結語
6.1本文從細觀力學角度，使用分塊映射網格劃分法結合各子塊網格隨機組合的思路，建立了混凝土三相細觀有限元模型，模擬了混凝土試件在單項拉伸時裂紋的產生、發展過程，結果與試驗結果定性相符。
6.2沒能模擬出應變軟化現象，綜合實驗、理論及數值方法對混凝土應變軟化現象的研究還需大量的工作。
參考文獻:
[1] 劉光廷，王宗敏. 用隨機骨料模型數值模擬混凝土材料的斷裂[J].清華大學學報，1996，36(1):84~89.
[2] 馬懷發，陳厚群，黎保琨.混凝土細觀力學研究進展及評述[J].中國水利水電科學研究院學報，2004，2(2) :124~130.
[3] 彭一江，黎保餛，劉斌.碾壓混凝土細觀結構力學性能的數值模擬仁[J ].水利學報，2001，6(10):19~22.
關鍵詞： 數值分析 課程教學 教學改革
“數值分析”又稱“數值計算方法”，是研究如何利用計算工具（如計算機）求解數學問題數值解的一門科學，因此數值分析隨著計算工具的發展而發展。它與傳統的數學課程有明顯的區別，它與概念、邏輯推理及各類問題求解的數學技巧相對而言不如先前數學課程那么重要，對具體問題求解方法的構造及利用數學理論對各類算法的收斂性和穩定性分析是課程的研究核心。根據這門課程的特點及授課對象的差異，本文提出了課程教學改革的思路。
一、更新觀念、改革課程內容
在理、工科數學類科目教學體系中，“數值分析”是一門重要的課程，理工科學生在各自專業學習中面臨各種專業問題需要建立數學模型求解，復雜的數學模型往往得不到解析解，那么如何求解？數值分析正好為他們提供了解決思路，因此能夠激發學生進一步學習數學、應用數學的意識和興趣，同時也能培養學生的創新思維和創新能力。
《數值分析》課程教學內容多，公式復雜，推導繁瑣，而學時少，教師不可能對每個知識點講解透徹，且算法實現對編程能力的要求較高，因此不少學生由于“掌握不了”容易厭學[1]。筆者認為根據學生不同專業、不同院校可以選擇符合自己特色的教材。該課程授課對象主要有兩類：一類是理科學生，他們學習的目標是不僅要會“使用”算法，還要會“創造”算法。另一類是工科學生，他們學習這門課程的目的主要是“使用”算法。由于這兩個專業群體學習的目標不同，教材的內容、體系及側重點應有差別。理科生學習時需要“研究”算法，學得相對精細點，學時不夠，因此在學習內容的選取上可以刪去一些在后續課程中將會學到的內容，比如常微分方程數值解。針對學生編程基礎薄弱，可選用基于MATLAB的數值分析教材[2]，由于MATLAB編程簡潔、數據處理方便，具有強大的數值計算功能，即使不具備較強編程能力的學生，學習這種教材也不至于為算法實現而發愁，為提高數值分析課程的教學質量創造了良好的條件。筆者認為基于MATLAB的數值分析教材符合創新性實踐教學的要求。當然對于編程能力強的部分理科學生來說，學習C語言等高級語言的教材有利于加深對算法思想的理解。
二、課程教學方法和教學手段的創新
學生是學習的主體，教師是教學活動的引導者、組織者，因此應充分發揮學生的主體作用。傳統的教學模式是教師在課堂上講授，學生被動接受，“滿堂灌”的教學模式，教師往往費力不討好。數值分析的教學目的不僅讓學生掌握算法的基本思想、基本方法和基本原理，還通過這門課程的學習提高學生數學素養，注重培養學生舉一反三的能力。因此，在教學內容的講解上要有所側重，以主帶次，在有限的學時里講清每一個主題，突出講授典型的、具有代表性并能體現其思想方法的常用算法和理論；對那些原理相近的內容不求面面俱到，可以精心設計教學提綱加以引導和提示，并且開發網絡教學平臺，把教學資源于網絡平臺，讓學生利用這一平臺以小組為單位課后討論自學完成。數值分析課程教學可以貫徹“少而精”的原則，達到創新教學的目的。
課堂教學應激發學生的興趣，在教學方式上，力求講清每一個主題的實際背景、目的和算法設計的出發點，通過對實際應用背景的描述，激發學生的學習欲望。不少學生在前期數學學習中感到困惑：數學到底有什么用？以實際問題為出發點，讓學生真切地體會到學習數值分析課程的意義所在。其次，充分調動學生的積極性和主動性，啟發學生進行獨立思考，帶領學生分析算法的適用范圍、優缺點，啟發學生對算法進行改進。讓學生明白任何一種算法都有局限性，有值得改進的地方，在課程學習中激發學生的科研興趣。針對某一個知識點，要求學生課前查閱資料認真預習，讓一些學生像老師備課一樣精心準備，并上臺講解。學生經過查閱資料、獨立思考這一主動學習過程提高興趣，感受到成功的喜悅，增強自信，由對這門課程最初的“害怕”逐步變為“喜愛”，也通過這一知識點的學習掌握學習方法。我們在信息與計算科學專業學生中進行這方面的嘗試，實踐證明通過這種以點帶面的學習，大大調動了他們的主動性和積極性。
三、重視實踐環節的教學和改革
傳統的課程教學重理論、輕實踐，而《數值分析》是一門實踐性很強的課程，學習這門課程時，應該將重點放在掌握數值分析的基本理論和思想及與結合計算機解決實際問題上，所以必須重視《數值分析》的實踐教學。
實踐教學環節可以分為兩大類：一類是課堂實驗，另一類是與工程實際背景相關聯的數學問題的求解，比如數學建模競賽、與學生專業領域相關的問題求解。對于課堂實驗，教師應編寫實驗指導書，鼓勵學生分組討論，讓學生通過實驗進一步體會算法思想，掌握算法精髓。比如學習用Newton迭代法解非線性方程時，要求初值比較接近真實解，這時引導學生實驗時選擇不同的初值進行比較，并歸納結論；由于實際問題是不知道真實解，難以滿足初值要求，因此對算法進行改進，提出了下山牛頓法，這一過程讓學生在實驗中自己感受，印象會深刻得多。對于工程實際問題的求解，可以引導學生分析全國大學生數學建模競賽歷屆試題，并將數值分析課程中的逼近思想，如插值法、最小二乘法、曲線擬合等方法用于求解中，我們必須將數學建模的思想融入數值分析的教學中。對于學生專業領域中的數學問題可引導學生分組討論，自己設計算法方案。通過實踐環節的教學，讓學生真正體會到學有所用，同時激發他們的學習熱情和科研興趣。
參考文獻：
關鍵詞：數值分析課程設計；實踐教學；教學改革；教學方法
中圖分類號：G420 文獻標志碼：A 文章編號：1002-2589（2013）24-0291-02
《數值分析》是一門研究數值求解各類數學問題的課程，是綜合性大學數學專業與某些理工科專業的必修課程，也是面向全校學生開設的通識選修課。在《數值分析》課程設計的實踐教學中，學生經常抱怨時間短，內容多且復雜，難學且不會應用。事實上，《數值分析》既有數學課程在理論上的抽象性和嚴謹性，又有解決實際問題的實用性和實踐性，它作為科學計算的基礎理論與基本方法，已經廣泛地應用于物理學、力學、計算機應用、航空航天、土木工程、機械工程、風險投資和經濟管理領域。如何改變學生的觀點、提高學生的興趣、提高《數值分析》課程設計的教學質量，已成為數學專業課程教育教學改革的焦點之一。筆者經過幾年的《數值分析》課程的教學，分析了《數值分析》課程設計實踐教學中存在的一些問題，并根據學生專業的特點提出幾點教學改革建議。
一、實踐課程中存在的問題
《數值分析》課程設計實踐教學，課程難學與本課程教學過程中的一些特點有關。具體說有如下幾點。
（一）任課教師的專業限制
《數值分析》課程是信息與計算科學專業的專業必修課，它的任課教師是數學專業教師。而《數值分析》課程設計培養的是厚基礎、寬口徑、應用型本科人才，屬于應用教育。在教學過程中，任課教師常常強調專業和基礎，忽視實踐與應用，這使得學生要花許多時間學習，才能學會該課程，如果不這樣做，很難達到老師的要求。這樣學生的實踐練習時間就很少，更難達到用所學知識來解決實際問題，造成教與學脫節、學與應用脫節。另外由于數學專業教師缺乏工科背景，教學中多是只給出概念、公式、定理，羅列出方法和計算，很難做到時時處處與實際應用需要結合，這樣會讓學生感覺該課程就是一門枯燥的數學課，很難做到學以致用，逐漸地喪失學習興趣。
（二）學生的專業特點限制
《數值分析》課程是面向數學系信息與計算科學專業開設的專業必修課，而信息與計算科學專業的學生數學基礎好，計算機基礎薄弱，程序操作能力差，課程設計實踐教學的目的是各種數值算法的計算機實現，學生本就覺得這門課程理論枯燥難學，由于不熟練編程序，很難有興趣自己操作實現各種算法，違背了課程設計的目的，使得這部分內容多是老師講解演練各種算法，學生很少動手操作，教學效果較差。
（三）數值分析課程設計實踐教學時間短，內容多
《數值分析》是一門與計算機密切結合的課程，該課程的任務是要根據計算機的特點給實際問題提供切實可行的算法，因而，課程設計是該課程教學中一個不可缺少的環節。然而，由于學時較少，在保證理論的限制下，課程設計實踐教學只有一周時間。而且《數值分析》課程設計多安排在教學的最后一周，學生在學期最后幾天往往會因很多事情影響學習，實際授課和實踐的時間就更少了，在這樣短的時間內完成許多算法的實踐是很難做到的，只能選擇幾個重點的算法由老師講解學生練習，學生只掌握了幾種算法的原理，掌握效果也并不好，很難用所學知識去解決實際問題，這與《數值分析》課程設計解決實際問題的實用性和實踐性相違背。
（四）數值分析課程設計教學模式單一
數值分析課程設計教學目標是教學生將所學的理論知識轉化為實踐能力，采用的教學模式就是上機講解和練習。連續五天的這種教學模式會讓學生覺得形式單調、內容枯燥、難學，在沒有老師在身邊監督下，有的學生干脆不聽不練習，玩起了游戲，老師過來檢查時關掉游戲，老師離開又繼續玩，反反復復地看著老師，幾天的實踐課程就這樣荒廢，很難達到預期的教學效果。
二、教學改革的幾點建議
針對上述的幾點問題，筆者在《數值分析》的教學中提出以下幾點改革建議。
（一）明確教學目的
必須明確我校人才培養的目標不是專業人才，而是應用型本科人才，培養操作能力強、寬口徑、復合型的人才。針對這種人才培養模式的特點，我們確立《數值分析》的教學目的是：通過這門課程的學習，使學生了解數值計算的基本原理，掌握一般的數值計算技能，增強學生實踐能力，使學生能夠在實際需要涉及無法求解的數學問題時，會使用數值計算的思想和方法，為以后的實際操作和研究工作打下良好的基礎。
（二）充實和精煉教學內容
《數值分析》是一門嚴謹、完整的學科體系，內容豐富。連同課程設計部分它的學時是64學時，筆者在近幾年的授課過程中感覺到，在少學時的情況下，如果所有內容都講到，往往會顧此失彼，不能突出重點，學生學下來后并沒有在腦海里對所學知識留下深刻印象。如果將實際應用較多的內容重點強調，相關的方法應用詳細講解，并在上機實驗上重點練習，其他知識點弱化講解，應用較少卻很繁瑣的內容精簡掉，這會使得與實際應用相關的重點內容有充足的時間講解，實用性很強的方法有動手操作的機會，不僅減輕了學生的學習負擔，還激發了學生思考用這些方法解決相關問題的興趣，增強了學生解決實際問題的能力。
（三）不斷改善實踐課教學模式和教學方法
《數值分析》研究的是各類數學問題的數值算法，對于實用性較強的算法一定要多給學生動手操作的機會，加深學生對這些算法的理解，熟練它們的應用。數值試驗的題目選擇一些與實際應用相關或生活相關的案例，讓學生將分析問題、建立模型、選擇合適的算法、編寫程序、分析結果這一計算主線完整地進行下來，以縮小學會該課程與解決實際問題的距離，增強學生的能力和信心。數學軟件可以指導學生使用簡單且方便調試的軟件，如Matlab等，這可以兼顧到學生計算機語言基礎薄弱，避免調試程序時語法錯誤頻出而喪失信心。在教學方法上，做以下改革：
1.采用在解決問題中引出理論內容的教學法
每節課前都提出一些與學生專業和生活聯系緊密的問題，和學生一起探討這些問題的解決方法，越簡單越好，從實際問題的解決方法入手，運用歸納、分析的手段引出相應的數值計算理論，之后從理論上研究解決問題的思想和方法，分析方法的優點、缺點以及所能解決問題的類型，進而給出解決實際問題的程序。例如在介紹數值積分時，課前不妨假設要購買一塊不規則形狀的養魚池，按面積付費，判斷所付價格是否合理。以其中兩點確定的方向為橫坐標軸，選擇坐標軸上的點后，用船行駛的時間乘以船速，得到不同點的縱坐標，每小塊面積用矩形面積近似，這些近似值加起來得到養魚池面積的近似，將這種思想和方法提煉出來就是數值積分的思想和方法。結合生活中的實例，學生既能掌握《數值分析》中各部分內容，又能了解其在生活中的應用，進而將遇到的問題分類，學會用不同的數值方法解決相應的問題。
2.采用優劣分析式教學方法
這里所涉及的優劣分析包含數值解法與解析解法優劣分析、同一問題不同解法優劣分析等。
如講到數值積分時，選幾個實際生活中遇到的問題為例，如購買性狀不規則的養魚池，需要估計面積，而在邊界曲線函數無法知道的情況下，不能直接應用計算定積分的牛頓-萊布尼茲公式來計算面積，而若被積函數是一張數據表、用牛頓-萊布尼茲公式積分后的原函數過于復雜等，都無法直接應用牛頓-萊布尼茲公式計算，而數值積分公式則可以克服這些弊端，強調了數值積分的實際應用價值，給出具體的數值算法，比較數值算法的計算量。又如，在介紹線性方程組的數值解法時，指出線性代數中的線性方程組的求解方法只對階數較低的方程組適用，階數稍高就會導致計算量很大，列舉一個五階線性方程組和一個十階線性方程組，觀察它們的計算量，而源自工程實踐中的線性方程組多是高階線性方程組，如土木工程中的有限元方法，最后都會歸結為解高階線性方程組，如果用傳統的解法會導致計算量相當大，而給出數值算法，比較數值算法與傳統算法的計算量，這樣，數值解法的優點可以吸引學生的學習興趣。
啟發學生思考身邊問題，發現問題、解決問題。課上，老師講授之外，師生共同圍繞學生提出的那些與生活相關的問題進行討論，提出解決問題的方法并在程序上實現問題的解決，對大家提出的算法進行點評，這種授課方式能夠激發學生學習的熱情。
（四）改革考核方法
以往的實踐課考核方式是給出幾道題目，給學生一定的時間，由他們自己設計算法，給出問題的解決方案，提交課程設計報告，根據這些報告給定成績。有的同學平時不練習，在最后抄襲別人的實踐報告。為了杜絕這類現象，在平時上課時，將學生分組，每次課給每組指定任務，大家分組討論，每組都給出問題的解決思路，每個組員都要給出補充，根據這些課上表現綜合給出考核成績，這既可以監督學生浪費課上時間玩游戲，又可以保證所給成績的公平性，調動學生積極性。
三、結束語
教學內容和教學方法的改革，是提高教學質量的關鍵，是高校提高自身的競爭力適應社會發展需求的必然趨勢。講授《數值分析》課程設計，必須做到與實際問題相結合，根據實際應用取舍教學內容，安排試驗學時，培養學生應用所學知識解決實際問題的能力，拓展學生思維，為培養基礎全面、操作能力強、綜合素質好的應用型本科人才奠定基礎。
參考文獻：
[1]李慶揚，王能超，易大義.數值分析：第5版[M].北京：清華大學出版社，2008.
[2]周生田，李維國.工科碩士研究生數值分析課程建設與教學改革[J].石油教育，2009，（1）：57-59.
關鍵詞：柴油機；螺旋進氣道；設計參數；缸內渦流比；數值分析
中圖分類號：TK422 文獻標識碼：A
Numerical Analysis of the Effects of the Design Parameters of a Diesel Helical Intake Port on the Swirl Ratio
HAN Zhiyu1,2，WANG Yong1，CHEN Zheng1，DENG Peng1
(1.Research Center for Advanced Powertrain Technology，Hunan Univ，Changsha，Hunan 410082，China; 2.State
Key Laboratory of Advanced Design and Manufacturing for Vehicle Body，Hunan Univ，Changsha，Hunan 410082，China)
Abstract:To study the effects of the design parameters of a helical intake port on the incylinder swirl ratio for an automotive diesel engine, the transit port and incylinder flows were simulated with the Converge CFD software. By changing the values of three main design parameters of the intake port, their effects on the incylinder swirl ratio were simulated and studied. The results have shown that the incylinder swirl ratio achieved the maximum values when the spiral chamber height, the vortex shell cutoff amount and the helical intake port angle were changed to the values of 13 mm, 1.27 mm and 0 degree, respectively. The results also indicated that the effects of the parameters on the swirl ratio were not additive when multiple parameters were varied at the same time. This study has shown that the numerical analysis of multiple design parameters of a helical intake port can be done with the use of the Converge software.
Key words:diesel engine；helical intake port；design parameters；incylinder swirl ratio；numerical analysis
隨著石油能源的不斷減少以及環境問題的日益加重，針對車用柴油發動機的燃燒和排放提出的要求也越來越嚴格.為了滿足日益嚴格的法規要求，必須合理的組織缸內燃燒過程，以提高燃料的燃燒效率和降低排放.
柴油機螺旋進氣道結構參數對缸內渦流的產生有至關重要的作用，合理的氣道結構參數不僅能保證在進氣過程中產生較強的渦流，還能保證在壓縮過程以及壓縮上止點附近缸內都具有較強的渦流.通過與噴油時刻的匹配，能保證在噴油燃燒時，缸內具有最優的渦流比，促進燃料與空氣的混合，提高缸內空氣的利用率，以達到改善發動機動力性、經濟性以及排放性的目的［1-2］.本文采用最新的三維CFD計算軟件Converge進行數值研究，分析了柴油機螺旋進氣道的3個主要結構參數對缸內渦流比的影響，實現了螺旋進氣道多結構參數的數值研究.
1 基本情況介紹
1.1 發動機基本參數
建立了一臺高速柴油發動機的氣道及燃燒三維模型，如圖1所示.氣道三維模型是通過三維激光坐標掃描儀掃描氣道砂芯獲取點云圖，并導入三維建模軟件中得到［3］.該柴油機基本參數如表1所示.
圖1 氣道及燃燒室三維模型
Fig.1 3D model of transit port and combustion chamber
表1 柴油發動機基本參數
Tab.1 Basic parameters of the diesel engine
型號
排量/L
壓縮比
燃燒室
形狀
進氣
方式
標定功率
/kW
最大轉矩
/(N?m)
直列
4缸
1.6
18.1
直噴
SymbolwA@
燃燒室
增壓
中冷
80(4 000
r/min)
230(2 000
r/min)
湖南大學學報(自然科學版)2012年
第4期韓志玉等：螺旋進氣道設計參數對渦流比影響的數值分析
1.2 三維模擬計算軟件介紹
本文采用了最先進的三維計算流體力學軟件Converge進行CAE研究.其不同于傳統CFD軟件之處，在于它的計算網格是在計算過程中自動生成，無需用戶提前畫出計算網格，從而為用戶節省了大量網格制作時間.在使用Converge時，用戶只需將制作的STL格式的幾何文件導入Converge的前處理軟件中進行簡單的幾何檢查和邊界設置即可［4］.
2 渦流比與進氣流量系數計算
2.1 渦流比與進氣流量系數
2.1.1 渦流比的計算
渦流比SR的計算公式為［4-5］：
SR=ΩflowΩcrankshaft. （1）
式中：Ωflow為缸內氣體的角速度，r/min；Ωcrankshaft為發動機的曲軸轉速，r/min.本文利用Converge計算軟件可直接計算出缸內渦流比的大小.
2.1.2 進氣流量系數
Ricardo無量綱流量系數C［6-7］F：
CF=QnAVV0,（2）
V0=2?ΔPρ, （3）
AV=π?d2V4. （4）
式中：Q為試驗測得的實際空氣體積流量，m3/s；n為進氣門數目；V0為理論進氣速度或速度頭，m/s；AV為氣門座內截面積，m2；dV為氣門座內徑，m；ΔP為進氣道壓差，Pa；ρ為氣門座處氣體的密度，kg/m3.本文通過Converge軟件得出計算的空氣體積流量Q，然后代入式（2），得出計算的流量系數CF.
3 邊界條件設定以及模型標定
3.1 邊界條件的設定
在發動機臺架上進行2000轉倒拖試驗并測得相關試驗數據.模擬發動機的工作范圍設定為進氣閥開啟前10°至壓縮上止點后20°，則對應的曲軸轉角范圍為-386°~20°.根據倒拖試驗測定的缸壓、溫度等，設定-386°時刻發動機的邊界條件如表2和表3所示.
表2 邊界條件1
Tab.2 Boundary condition 1
進氣溫度
/K
進氣壓力
/kPa
缸內溫度
/K
缸內壓力
/kPa
排氣溫度
/K
排氣壓力
/kPa
325.8
124.2
352
245.2
342.8
196.5
表3 邊界條件2
Tab.3 Boundary condition 2
活塞溫
度/K
氣缸溫
度/K
缸蓋溫
度/K
進氣閥
溫度/K
進氣道
壁面溫
度/K
排氣閥
溫度/K
排氣道
壁面溫
度/K
380
353
353
353
343
370
370
3.2 模型的標定
將數值計算的發動機缸壓與倒拖試驗所得到的實驗缸壓進行對比，如圖2所示.
由圖2可以看出：在開始時刻計算缸壓和實驗缸壓存在微小的波動，在之后的時刻，計算缸壓和實驗缸壓十分吻合.在壓縮上止點時刻，計算缸壓略高于實驗缸壓，計算的最大缸壓（5.812 1 MPa）比實驗的最大缸壓（5.718 4 MPa）高1.6%.因此，可以認為：數值模型能很好地反映實際發動機的工作情況，可以使用此數值模型進行下一步研究.
曲軸轉角/° ATDC圖2 計算缸壓與試驗缸壓對比
Fig.2 Comparison between calculation and
experimented results
4 調整結構參數的數值模擬研究
4.1 結構參數介紹
將螺旋進氣道的三維模型在CAD軟件中直接導出二維工程圖并進行適當簡化，得出決定螺旋進氣道結構形狀的主要結構參數，如圖3所示.
圖3 螺旋進氣道主要結構參數
Fig.3 Main structure parameters
of the helical intake port
結構參數包括：渦殼切割量η、氣道入口截面積S1、氣道最小截面積S2、螺旋底坡角β1、螺旋坡角β2（氣道最小截面的法線與水平面的夾角）、螺旋室高度H、螺旋室直徑D、螺旋進氣道高度μ、渦殼半徑R1、R2、螺旋進氣道轉角θ、螺旋進氣道偏心距e以及渦殼轉角β3.本文主要研究螺旋室高度H、渦殼切割量η以及螺旋進氣道轉角θ改變對缸內渦流比的影響.
4.2 參數H的影響
如圖4所示，H為螺旋進氣道的螺旋室高度，通過三維造型軟件改變H的取值，分別取H1=7 mm，H2=7.5 mm，H3=10.5 mm，H4=13 mm以及H5=13.5 mm.計算得出H不同取值時的缸內渦流比，如圖5所示.對比H不同取值的進氣量，如圖6所示.由圖5可以看出，在整個壓縮過程中，缸內渦流比呈先減小后在上止點附近快速增大，達到一個峰值以后再減小的過程［8］.原因是：在壓縮過程中，隨著活塞的上移，進氣初期形成的缸內渦流受壓而減弱，從而使渦流比降低；當接近上止點時，缸內氣體被強制壓入燃燒室，氣體旋轉半徑減小并產生擠流，從而使缸內渦流增強，渦流比增大并在上止點附近達到最大值；隨著活塞下移，缸內氣體從燃燒室內流出，氣體旋轉半徑再次增大，缸內渦流減弱，渦流比減小.
圖4 螺旋進氣道結構參數H
Fig.4 Helical intake port structure parameter H
曲軸轉角/°ATDC圖5 H不同取值時的渦流比
Fig.5 Swirl ratio of different values of parameter H
曲軸轉角/°ATDC圖6 H不同取值的進氣量
Fig.6 Air inlet of different values of parameter H
同時，可以看出隨著H取值的增加，缸內渦流比整體上成先減小后增大再減小的規律，并在H=13 mm時，缸內渦流比獲得最大值.可以得出：H與缸內渦流比之間存在非線性的相關關系，存在一個最佳H取值，使渦流比具有最大值.
繼續分析圖5可知，改變H的取值，對-20°和0°時刻的渦流比影響較大，而對20°時刻缸內渦流比的影響不明顯.因此，-20°和0°時刻的缸內渦流比對結構參數的變化較敏感，結構參數的微小變動，都會影響此時的缸內渦流比，而20°時刻的缸內渦流經過了上止點的壓縮過程，對結構參數的變化已經不那么敏感了.
由圖6可知，雖然H的取值不斷增加，但進氣量的變化卻不明顯.在-120°時刻，缸內進氣量最大值與最小值之間相差0.72%，因此，可以認為：結構參數H的改變對進氣量沒有影響.
4.3 參數η的影響
根據實際試驗中，渦殼切割量對缸內渦流比有一定的影響，因此，取η為渦殼切割量，如圖7所示.通過修改使其取值分別為η0=0 mm，η1=1.27 mm，η2=2.54 mm.
圖7 結構參數η
Fig.7 Structure parameter η
計算得出η不同取值時的缸內渦流比，如圖8所示.對比η不同取值的進氣量，如圖9所示.
曲軸轉角/°ATDC圖8 η不同取值的渦流比
Fig.8 Swirl ratio of different values of parameter η
由圖8可以發現，隨著η的增大，缸內渦流比先增大后減小.在η=1.27 mm時，渦流比在-180°~20°曲軸轉角之間都具有最大值，且其最大渦流比（1.114 8，-4°）比η=0 mm的最大渦流比（0.950 6，-3°）高17.3%；同時，在-20°和0°時刻，η=1.27 mm的渦流比比η=0 mm時高17.9%和17%.隨著η值的繼續增大，缸內渦流比降低，并低于η=0 mm的渦流比.
曲軸轉角/°ATDC圖9 η不同取值的進氣量
Fig.9 Air inlet of different values of parameter η
由圖9可知，雖然η的取值不斷增加，但進氣量的變化卻不明顯.在-120°時刻，缸內進氣量最大值與最小值之間相差0.58%，因此，可以認為：結構參數η的改變對進氣量沒有影響.
4.4 參數θ的影響
以進氣閥的軸線為軸，進氣閥的軸心與進氣道入口截面的中心連線，繞此軸旋轉的角度為θ，如圖10所示.設氣道初始位置的θ角度為0°，且順時針旋轉為+，逆時針旋轉為-.通過三維建模軟件修改θ的取值，分別為-15°，-10°，-5°，+5°，+10°以及+15°.
計算得出缸內渦流比，如圖11所示.對比參數θ不同取值的進氣量，如圖12所示.由圖11可以看出，無論是順時針旋轉還是逆時針旋轉，缸內渦流都降低，且順時針旋轉降低的程度要遠小于逆時針旋轉降低的程度，同時隨著旋轉角度的增大，渦流比下降的速率也逐漸增大.根據不同的θ取值，可以得出初始位置時螺旋進氣道所處的θ角是使缸內渦流比為最大值的最佳角度.
圖10 結構參數θ
Fig.10 Structure parameter θ
曲軸轉角/°ATDC曲軸轉角/°ATDC
圖11 θ不同取值的渦流比
Fig.11 Swirl ratio of different values of parameter θ
曲軸轉角/°ATDC圖12 θ不同取值的進氣量
Fig.12 Air inlet of different values of parameter θ
從圖12可以看出，無論θ的取值如何變化，進氣量的變化卻不明顯.在-120°時刻，缸內進氣量最大值與最小值相差0.6%，因此，可以認為：結構參數θ的改變對進氣量沒有影響.
5 調整多個結構參數的研究
為研究多個參數改變時，對缸內渦流比的影響，分別同時改變2個參數和3個參數，設置5個Case的模擬計算，表4所示，為每個Case中各參數的取值.
表4 各Case參數取值
Tab.4 The parameters values of different Cases
名稱
H/mm
η/mm
θ/（°）
Case1
10.5
Case2
13
1.27
Case3
12
1.27
Case4
13
1.27
+10
Case5
13
1.27
-10
其中，Case1的參數取值為所研究發動機螺旋進氣道結構的實測值；Case2的參數取值是將各參數調節到使渦流比為最大值的參數取值；Case3的參數取值是保持Case2中η和θ的取值不變，調整H的取值；Case4和Case5的參數取值是保持Case2中H和η的取值不變，調整θ的取值.
按表4所示，修改結構參數后導入計算，得到5個Case的渦流比，如圖13所示.對比5個Case的進氣量，如圖14所示.
曲軸轉角/°ATDC圖13 5個Case的渦流比
Fig.13 Swirl ratio of five Cases
曲軸轉角/°ATDC圖14 5個Case的進氣量
Fig.14 Air inlet of five Cases
比較圖14中的Case1和Case2可以發現，雖然分別單獨調節H=13 mm或η=1.27 mm時能使渦流比達到最大值，但當同時調整H和η時，并沒有使渦流比得到提高.因此，可以認為：各結構參數對渦流比的影響效果不是簡單的相互疊加，而是存在相互制約、相互抵消的作用.在此基礎上，對比分析Case2和Case3，在保持η和θ取值不變的情況下，適當調整H的取值時，發現Case3的渦流比得到一定幅度的提高，但整體上仍小于Case1的渦流比.同樣，對比分析，Case2和Case4和Case5可以發現，在保持H和η取值不變的情況下，調整θ的取值，渦流比較Case2有增大也有減小，但變化幅度沒有調整參數H明顯.
通過上面的分析，可以得出：缸內渦流比是各結構參數相互作用的結果，各結構參數之間存在相互制約的作用.欲使缸內渦流比得到優化，必須準確了解各結構參數與渦流比的相關關系，通過分析相關關系可以得出優化缸內渦流比的各參數取值.
分析圖14可以看出，在各種情況下進氣量的變化都不大.在-120°時刻，缸內進氣量最大值與最小值相差1.1%，因此，可以認為：在誤差允許的范圍內，本文中所討論的結構參數改變對進氣量沒有影響.
6 缸內流動情況的分析
分析修改結構參數后對缸內流動情況的影響，沿兩進氣閥的軸線做一豎直截面，取進氣閥在最大升程時刻，且右側為螺旋進氣道，獲得缸內流動速度分布圖，如圖15所示.
(a)Case1(b)Case 5圖15 缸內流動速度分布
Fig.15 Distribution of incylinder flows speed
分析圖15可以看出：在15(a)圖中，氣缸中心處產生了一個較大渦流，并進一步使進氣氣流在兩氣門之間向左側流動，加速了缸內渦流的產生；而15(b)圖中，由于缸內產生的渦流較小，進氣氣流進入缸內后直接流向燃燒室，缸內不能很好地組織渦流.分析原因，可以得出：由于Case1和Case5的螺旋進氣道在結構設計參數H，η和θ的取值上存在一定的差異，從而使Case1和Case5的進氣道在進氣過程中產生渦流的能力發生改變.Case1中各結構設計參數的取值優于Case5，從而使Case1的缸內渦流比高于Case5.因此，合理的設計螺旋進氣道的結構參數能很好的改善發動機的缸內流動情況，提高缸內渦流比，最終達到優化發動機缸內燃燒過程的目的.
7 結 論
1）本文討論的螺旋進氣道結構參數的改變對進氣量的影響較小，而對缸內渦流比的影響較大.當分別單獨修改螺旋室高度、渦殼切割量以及螺旋進氣道轉角為13 mm，1.27 mm以及0°時，缸內渦流比得到大幅度的提高，并取得最大值.
2）分析了發動機缸內渦流比的動態變化情況，研究了對燃燒有意義的時間段（-20°，0°和20°）的缸內渦流比，其結果比穩流試驗臺的結果更反映發動機的實際工作過程，更真實合理.
3）同時改變多個結構參數對缸內渦流比的影響比較復雜，各參數之間的作用效果不是簡單的相互疊加的，而是相互制約的.
4）采用三維計算流體力學軟件Converge，可以較方便地實現螺旋進氣道多結構參數的數值研究和優化，對氣道的設計有重要的指導意義.
致謝
感謝Convergent Science, Inc.提供Converge軟件以及相應的技術支持.
參考文獻
［1］ 周龍保．內燃機學［M］．3版．北京：機械工業出版社，2011：261-266．
ZHOU Longbao．Internal combustion engine theory［M］．3nd ed．Beijing：China Machine Press，2011：261-266．(In Chinese)
［2］ 蔣德明．內燃機燃燒與排放學［M］．西安：西安交通大學出版社，2001：512-538．
JIANG Deming．Internal combustion engine combustion and emissions［M］．Xi`an：Xi`an Jiao Tong University Press，2001：512-538.(In Chinese)
［3］ 王天友，劉大名，沈捷，等．四氣門柴油機進氣道開發［J］．內燃機工程，2008，29（2）：51-55．
WANG Tanyou，LIU Daming，SHEN Jie，et al．Development of Intake ports of foulvalve diesel engine［J］．Chinese Internal Combustion Engine Engineering，2008，29（2）：51-55.(In Chinese)
［4］ RICHARDS K J，SENECAL PK, POMRANING E.CONVERGE (Version 1.3)［M］.Middleton：ConvergentScience，Inc,2008．
［5］ 吳克啟，舒朝暉．高等流體力學［M］．北京：中國電力出版社，2009：13-26．
WU Keqi，SHU Chaohui．Advanced fluid mechanics［M］．Beijing：China Electric Power Press，2009：13-26．(In Chinese)
［6］ 周俊杰，徐國權．FLUENT工程技術與實例分析［M］．北京：中國水利水電出版社，2010：14-17．
ZHOU Junjie，XU Guoquan．FLUENT engineering techniques and examples analysis［J］．Beijing：China WaterPower Press，2010：14-17(In Chinese)．
［7］ 邱卓丹，沈捷．直噴式柴油機螺旋進氣道性能試驗及評價方法［J］．內燃機工程，2005，26（3）：27-30．
摘 要:根據現場地質條件，以地基關鍵斷面為研究對象建立相應的數值模型，應用自行開發RFPA2D程序就廣州地鐵二號線盾構開挖對新南方購物中心大樓穩定性影響問題進行了數值模擬研究·分析了隧道斷面洞頂的沉降位移和建筑物地基的沉降位移，數值模擬結果再現了隧洞開挖引起的周圍巖體的損傷破壞過程，得到了隧洞及地基周圍巖土的應力變化圖，塑性區破壞圖，以及位移矢量圖，預報了破碎帶的未來發展趨勢·數值計算結果與采用解析算法的結果進行了對比驗證，結果與現場觀測資料基本吻合·
關鍵詞:地鐵開挖;地基沉降;數值模擬;建筑物穩定性;有限元法;塑性區;巖體損傷
由于城市地鐵開挖引起的相應的問題(比如說地基沉降、地表沉陷等)近年來在實際工程中顯得尤為關注，特別是地鐵開挖引起地基沉降的位移變化情況、應力分布狀態以及隧洞的穩定性是設計者和施工單位所關心的·為了限制巖土體的變形破壞，必須選用適當的開挖方法和相應的加固措施·正確預測開挖后地鐵周圍巖土體的位移變化、應力矢量分布以及隧洞的變形破壞和破壞范圍顯得十分重要·但由于工程開挖問題涉及到的地質因素、工程因素復雜[1，2]，本文提出采用自行開發的RFPA2D數值模擬技術和工程分析法對隧洞盾構施工過程引起淺基礎下土層的應力應變進行分析·
工程實例以廣州地鐵盾構開挖對新南方購物中心大樓穩定性影響的分析為例·
1 工程概況
廣州地鐵二號線客村—鷺江站區間是盾構法的施工區間·該區間隧洞穿越六層的新南方購物中心時，有一排柱的基礎位于開挖隧洞的側上方，垂直距離是7.79~8.69m·按設計單位提供的地質資料、地鐵斷面寬度為10m，最高拱高為6m，平均埋深為11m，圍巖主要為砂巖，具體地質及隧洞相關位置如圖1所示，地層由上至下分布順次為:雜填土(1)、粉質粘土(2)、硬塑，標貫擊數N為16~36，強風化泥質砂巖(3)、中風化泥質砂巖(4)、微風化泥質砂巖(5)的層狀分布·
關鍵詞：拱形支護結構；空間特性；內力；變形
Based on the MATLAB numerical analysis of arched retaining structure is analyzed
Li Xin
（HeBei Ruizhi Transportation Technology Consulting Co., Ltd
HeBei Provincial Engineering Eesearch Center for Road Structure and Material
Shijiazhuang050091）
Abstrct: Arched retaining structure of deep foundation pit can give full play to its space effect, the uniform stress, the cross section only compressive stress, which can effectively control the deformation and ensure the stability of foundation pit. Engineering on the arch supporting system, using the shell theory to establish the differential equations, and then calculate the internal force and deformation of retaining structure, due to considering the spatial characteristics of this kind of supporting structure, so this method is more in line with the actual working situation.
Key words: arch supporting structure; spatial features; the internal force; deformation
中圖分類號： TV551 文獻標識碼： A
1前言
拱結構擁有良好的力學性能，可以將上部荷載沿著拱圈傳遞至拱腳，從而在拱圈內拉應力值很小，甚至不出現拉應力[1~2]。拱形支護結構在基坑開挖支護中應用較為廣泛，這種結構體系與常規的直線形支護結構相比，具有較多優點，工程實踐已表明采用這種支護形式后可大大減低基坑支護的造價，經濟效益相當可觀[4~5]。
2結構計算模式
拱形支護結構通常由幾段曲率不同的拱組成，在每個拱與其他結構連接處（拱腳）受力最大，故拱腳處往往得到加強，每一段拱從空間上看都可認為是直線邊支承的殼體結構，它所受的荷載為土壓力，方向為殼體的法向，見圖1所示。圓柱殼的內力分布形式見圖2所示。在基坑工程中，圓拱使用較多，現以圓柱殼作為支護結構的計算模型。對于只承受法向荷載的圓柱殼，其基本微分方程為式（1）：
（1）
其中：分別為圓柱殼的母線、環向坐標；分別為圓柱殼在母線、環向及法向上的位移；E為殼體材料的楊氏模量；為殼體材料的泊松比；R為圓柱殼的半徑；t為圓柱殼厚度；P為殼體所受的法向荷載（即土壓力）。
圖l拱形支護結構計算模型
圖2（a） 圖2（b）
圖2拱形支護結構計算模型
引入位移函數F（），則可將殼體的位移表示成式（2）：
（2）
將式（2）代人式（l），前兩式滿足，由第三式可得：
（3）
式中，
對于拱形支護結構，把拱腳處假設為簡支端，則位移函數F（）可表示為單三角級數形式：
（4）
式中，m為大于1的整數。
將法向荷載P（）也展開成三角級數形式為
（5）
將式（5）及式（4）代人式（3），比較方程兩邊的系數，可得
（6）
式中，。
對于拱形支護結構，被動區的土壓力可以用彈簧模擬，即，其中為水平基床系數；w為支護結構的位移。
主動區的土壓力采用主動土壓力，因此對于拱形支護結構其主動土壓力在母線方向上是線性變化的，而在環向上保持不變。
殼體支護結構上作用的荷載為：P=，其中h為基坑深度；c，d為與土壓力有關的量。
支護結構在坑底以上部分的微分方程為：
（7）
殼體的內力由位移函數表達為：
根據基坑開挖面處（a=h）支護結構的連續性，有：
位移連續：
力的連續：
其中上標為u的符號表示開挖面以上的量，上標為d的符號表示開挖面以下的量。根據基坑頂部和底部的邊界條件，得：
邊界為固支時：；
邊界為間支時：；
邊界為自由端時：，。
由連續條件及邊界條件即可求出式（4）的位移函數，根據式（2）和位移函數表達式可計算出支護結構各點的位移及內力。
3工程實例
某深基坑為局部扇形結構，開挖深度6.8m，采用由單排650鉆孔灌注樁組成的拱形支護結構，圓拱半徑為25m，圓心角為120o，在兩邊拱腳處各布置10根650鉆孔灌注樁，樁側加打水泥攪拌樁以加固周圍土體，形成支撐體系。支護結構平面布置及剖面見圖3所示。該深基坑主要土層自上而下分布情況為：雜填土，厚度1.5m；粘土，厚度1.0m；淤質粘土，厚度11.0m；粉質粘土，厚度3.2m；黃褐色粘土，厚度4.5m。為了監測基坑附近土體的側向位移，在拱頂處設置測斜管進行現場實時實測。
圖3支護結構平面及剖面
根據土層分布及地下水位情況，計算得到土水壓力分布情況。利用MATLAB數值分析，得到拱頂處的位移（w）計算值與實測結果，見圖4所示；拱頂彎矩（M）計算結果，見圖5所示。監測結果表明，盡管為單排樁組成的自立式支護結構，但由于采用了拱的結構形式，使整個支護結構的側向位移和內力均較小。
圖4拱頂處側向位移計算域實測值比較
圖5拱頂處彎矩計算值
由圖3可知，拱頂處的側向位移較小，最大值僅為3.6cm，側向位移在0~6m范圍內隨深度的增加而增加，6m處的側向位移最大；6m以下范圍內側向位移隨深度的增加而減小。由圖4~圖5可知，拱頂處的彎矩絕對值也較小，最大值僅120kN.m，0~5m范圍內負彎矩絕對值隨深度增加而增大，5m以下范圍負彎矩絕對值隨深度增加而減小，10~15m之間出現反彎點，反彎點以下范圍內正彎矩隨深度增加而增大。
4結語
通過上述分析，得到以下幾個方面的結論與建議：
（1）拱形支護結構是一種合理的支護形式，在本文介紹的非閉合情況下，拱頂處的側向位移及內力均較小，同時又充分發揮材料的抗壓性能，比常規直線形支護結構具有較大的優越性。
（2）拱形支護結構可以把支護形式由單拱向多拱方向發展，推動拱形支護結構向更大、更高結構層次發展和應用。在材料選用上，可以向更為廣泛的鋼筋混凝土灌注樁、地下連續墻、鋼樁支護等支護方式上拓展拱形結構，有效降低工程成本、節約工期。
（3）由于拱形支護結構固有的外形特點，使得拱形支護結構在施工場地狹小的工程在使用上受到限制。
參考文獻：
[1] 陸瑞明,趙錫宏.擋土拱圈的內力和變形計算[J].巖土工程學報,1997,19(6).
[2] 張祖聞,陳新余.高層建筑地下結構及基坑支護[M].宇航出版社,1994.
[3] 吳連元.板殼理論[M].上海交通大學出版社.1989,12~15
[4] 梁春江.連拱式基坑支護與實踐[J].施工技術.2002.
[5] 房集明.連拱形基坑支護體系[J].巖土工程技術.2000
[6] 龔曉南等.深基坑工程設計施工手冊[M].中國建筑工業出版社.1998
[7] 張士喬等.基于拱梁法原理的深基坑拱形圍護結構分析[J].土木工程學報.2001
[8] 黃生根.連拱式水泥土支護結構的應用[J].巖土工程技術.1998
作者詳細信息：李欣，女，（1986―），助理工程師；
