Monte_Carlo方法必備知識

1、目錄O Monte Carlo方法簡介d Potts模型Monte Carlo方法的應(yīng)用。課堂思考題材料設(shè)計優(yōu)化與生物醫(yī)用材料研究室Timn Length(year)io7(Oebye)1O°large scale FEN4, FD<largo-scale plasticity elasticity heat t ran sporty etc.)cellular automata, percolation models (scale indepenctent)crystal plasticity FBMO O a ohomogeniztion methods: Taylor,

2、Sachs. Voigt, Reuss, Hashin-Shtrikman, and self- consistent models in conjunction with advanced constitutive lawsdynamicsG i n zburg-Land au-ty pe and kinetic phase field models, microscopic field kinetic models(mm<T multi state Potts modelstopological network and vertex models, boundary dynamics

3、spring modelsmolecular dynamics(E)Metropolis Morrte Carlolocal electron density functional theory(Raabe D,計算材料學(xué),北京,化學(xué)工業(yè)出版社，2002 )Figure 1.2: Typical space and time scales and simulation methods in computational materials science 模擬計算方法介觀至宏觀尺度:有限元法、邊界元法有限差分法、離散單元法微觀至介觀尺度：Potts模型Monte Carlo方法相場動力學(xué)方法(P

4、hase-field model) 兀胞自動機方法(Cell Automation) 表面演化法(Surface evolver) 頂點法(Vertex model)納觀至微觀尺度：分子動力學(xué)方法Monte Carlo 方法一 Monte Carlo方法簡介 L蒙特卡羅(Monte Carlo )方法這一方法作為一種模擬工具，其系統(tǒng)發(fā)展始于1944年。馮諾伊曼 (J.von.Neumann)> N. Metropolis.烏拉姆(S.Ulam)在其參加的"曼哈頓計 劃”研究時引入了蒙特卡羅方法，他們用其模擬了中子在可裂變材料中 的無相關(guān)空間擴散(即無規(guī)行走問題)。同時邁耶(J.

5、E.Mayer)也獨立的 提出了一種類似的數(shù)值方法，通過非權(quán)重隨機抽樣過程用其計算積分問 題。人們認識到確定性問題可以通過恰當?shù)母怕誓Ｐ蛠砗苋菀椎慕鉀Q。Monte Carlo模擬通常是采用無相關(guān)隨機數(shù)進行大量的計算機實驗。 鑒于大量運用隨機抽樣，Monte Carlo方法的發(fā)展是與計算機技術(shù)的進 步密切相關(guān)的。(Raabe D,計算材料學(xué)，北京，化學(xué)工業(yè)出版社，2002 )蒙特卡羅方法概述：一種隨機模擬(random simulation)方法,有時也稱作隨 機抽樣(random sampling)技術(shù)或統(tǒng)計實驗(statistical testing)方法。以概率和統(tǒng)計理論方法為基礎(chǔ)的一種計

6、算 方法，是使用隨機數(shù)(或更常見的偽隨機數(shù))來解決很多 計算問題的方法。將所求解的問題同一定的概率模型相聯(lián) 系，用電子計算機實現(xiàn)統(tǒng)計模擬或抽樣，以獲得問題的近 似解。1777年比豐針問題取一根長度為1(1 <d)的針,隨機地向畫有間距為d的多條平行 線的平面上擲”次，記針與直線相交的次數(shù)為加。經(jīng)計算可 以知道，針與直線相交的概率約為：m 21 實驗者年份投計次數(shù)沫爾彈 wol巧I 850500()Wr 呻f JSXS m ith)1S5532043.1 5531 W9斗1 U 203/14 W拉査里尼 <Lazarini)1-90134083.141592從表中數(shù)據(jù)可以看到，一直到

7、公元20世紀初期，盡管實驗次數(shù)數(shù)以千 計，利用蒙特卡羅方法所得到的圓周率，還是達不到公元5世紀祖沖 之的推算精度。這可能是傳統(tǒng)蒙特卡羅方法長期得不到推廣的主要原 因。計算機技術(shù)的發(fā)展，使得蒙特卡羅方法得到快速的普及。現(xiàn)代的蒙特 卡羅方法，已經(jīng)不必親自動手做實驗，而是借助計算機的高速運轉(zhuǎn)能 另，使得原本費時費力韻實驗過程，變成了侯速和輕而易舉的事情。 它不但用于解決許多復(fù)雜的科學(xué)方面的問題基本思想：當所求解問題是某種隨機事件出現(xiàn)的概率，或 者是某個隨機變量的期望值時，通過某種“實驗” 的方法，以這種事件出現(xiàn)的頻率估計這一隨機事 件的概率，或者得到這個隨機變量的某些數(shù)字特 征，并將其作為問題的解。

8、一般步驟（1）構(gòu)造或描述概率過程對于本身就具有隨機性質(zhì)的問題，如粒子輸運問題，主要 是正確描述和模擬這個概率過程，對于本來不是隨機性質(zhì)的確 定性問題，比如計算定積分，就必須事先構(gòu)造一個人為的概率 過程，它的某些參量正好是所要求問題的解。即要將不具有隨 機性質(zhì)的問題轉(zhuǎn)化為隨機性質(zhì)的問題。（2）實現(xiàn)從已知概率分布抽樣構(gòu)造了概率模型以后'由于各科概率模型都可以看作是由 各種各樣的概率分布構(gòu)及的，因此盧王已跖概車分布的隨機變 量（或隨機向量），就成為實現(xiàn)蒙特卡羅方法模擬實驗的基本 手段，這也是蒙特卡羅方法被稱為隨機抽樣的原因。最簡單、 最墓本、最重奏胡一個彳既率分希是（d 1 ）上的均勻分希（

9、或 稱矩形分布）。產(chǎn)生隨機數(shù)的問題'就是從這個分布的抽樣問 題。隨機數(shù)是我們實現(xiàn)蒙特卡羅模擬的基本工具。(3 )建立各種估計量一般說來，構(gòu)造了概率模型并能從中抽樣后，即實現(xiàn)模擬 實驗后，我們就要確定一個隨機變量，作為所要求的問 題的解，我們稱它為無偏估計。建立各種估計量，相當 于對模擬實驗的結(jié)果進行考察和登記，從中得到問題的 解。蒙特卡羅方法在金融工程學(xué)、宏觀經(jīng)濟學(xué)，生物醫(yī) 學(xué)，計算物理學(xué)（如粒子輸運計算、量子熱力學(xué)計 算、空氣動力學(xué)計算）等領(lǐng)域應(yīng)用廣泛。 Monte Carlo方法在材料學(xué)中的應(yīng)用主要涉及到表面 與界面、擴散、相變、聚合物、外延生長及結(jié)構(gòu)、異 相界面、晶體生長、斷裂及

10、材料熱力學(xué)性質(zhì)等方面。二Potts模型Monte Carlo方法的應(yīng)用Potts模型Monte Carlo方法在材料微結(jié)構(gòu)模擬中具有廣泛 的適用性，并且已經(jīng)應(yīng)用于靜態(tài)初次再結(jié)晶、正常晶粒長 大、二次再結(jié)晶、動態(tài)再結(jié)晶、有粒子存在的再結(jié)晶、凝 固和擴散等領(lǐng)域。我們MMMD課題組也對Potts模型 Monte Carlo方法的改進和應(yīng)用進行了 一定的研究。下面結(jié)合本組的工作，主要以單相材料正常晶粒長大為例，對Potts模型Monte Carlo方法進行介紹。Gain Boundary Dynamics as a Tool for Microstructure ContolMaterial Prop

11、ertiesMotion晶粒長大：無應(yīng)變多晶體材料在退火過程中系統(tǒng)平均晶粒尺寸逐漸增 大的現(xiàn)象。晶粒長大可以是初次再結(jié)晶的后繼過程，即發(fā)生于形變試 樣初次畀結(jié)晶完成以后的繼續(xù)退火過程中，也可以發(fā)生在無鷹始形變 試樣的退火處理過程中。晶粒長大可以分為正常晶粒長大和異常晶粒 長大。正常晶粒長大的特點是長大速度比較均勻，在長大過程中晶粒的尺寸 分布和形狀分布幾乎不變。異常晶粒長大是組織中少數(shù)晶粒吞并基體 中其他較小的晶粒而長大。某種意義上講，晶粒長大研究是一個金屬學(xué)理論問題，但就這一研究 的起源和最終服務(wù)目的而言，晶粒長大研究是與材料性能密切相關(guān) 的。隨著人們對材料的組織、結(jié)構(gòu)與性能之間相互關(guān)系認識

12、的深入， 越來越顯出晶粒長大研究對控制和改善材料性能的重要性。研究晶粒長大的目的之一是控制晶粒尺寸。晶粒尺寸既反 映金屬材科的微觀組織特征，又直接影響材料的性能。例 如低碳鋼中晶粒尺寸與才扌料的機械性能、脆性轉(zhuǎn)變溫度有 直接關(guān)系。1. 細化晶粒結(jié)構(gòu)鋼：改善韌性同時提高強度變形鋁合釦提高強度，改善產(chǎn)品表面粗糙度和提高變形能力超塑性合金：提高其常溫強度而降低其高溫強度，實現(xiàn)超塑性的關(guān)鍵2 避免晶粒過細高溫合金、某些鎂鋁合金（防止產(chǎn)生應(yīng)變痕跡）等3.控制材料的織構(gòu)組態(tài)深沖加工的板材：為防止產(chǎn)生制耳，需防止過分強烈的單一織構(gòu)硅鋼片：為獲得高的磁導(dǎo)率，降低能耗，需控制織構(gòu)組態(tài)3D晶粒長大動畫MCS=10

13、05000, step=100作者:由于不透明材料3D微觀組織的不直接可視性，導(dǎo)致許多重要的 3D組織表征參量無法直接測量，許多涉及3D顯微組織演變的材 料理論模型無法得到驗證O采用Potts模型Monte Carlo方法研究晶粒長大，則可以有效避 免上述缺點。利用這種既遵從材料顯微組織形成和演變規(guī)律， 又已數(shù)字化且可視化的顯微組織仿真的靜態(tài)或動態(tài)模型，可以 進行晶?；蛉魏谓M織組成物及其動態(tài)演變過程的直觀分析和定 量研究（Exner教授將其稱為“虛擬金扇學(xué)T，獲得若干真實金 相學(xué)/斤無法獲得的組織表征彳言息和含時間變量陽動力學(xué)顯微組 織數(shù)據(jù)，將有助于我們對真實材料顯微組織及其各種演變過程 的進

14、一步了解。Exner H E. Practical Metallography,1999;36(3):115-137材料設(shè)計優(yōu)化與生物醫(yī)用材料研究室材料設(shè)計優(yōu)化與生物醫(yī)用材料研究室顯微組織形態(tài)a圖1 Radhaki-islman改進Monte Carlo算法晶粒長大過程不同階段的兩個多晶體組織三維組織模型晶粒數(shù)目分別為a)10085和b)2009(Qin X. PhD Dissertation. USTB, Beijing2003? p.l 13)1、標準Potts Monte Carlo模擬的原理及方法晶粒組織的晶界在本質(zhì)上是取向不同的晶粒之間的界面，而晶粒又是 處于晶格點陣位置上的原子所組

15、成的。如果將原子所處晶格點陣的取 向（即原子所屬晶粒的取向）稱為該原子的“取向?qū)傩?則晶粒組織 中的晶界可理解為是由取向?qū)傩圆煌脑又g的交界點連接構(gòu)成555 55 5 55 51111111111111111111111111111(a)222222222 2 2 2 22 2 2 22 2 2 21 1 1 12 2 21 1 1 2 2 2 皿2 2 2 21 1 1 q1 2 2 2 2 222222222波茨模型中晶界表征示意圖退火過程中晶界的移動即晶粒長大過程可理解為晶界附近原子按一 定規(guī)律改變其取向?qū)傩运鸬摹昂暧^”效應(yīng)。假設(shè)模擬系統(tǒng)中有N個原子（點陣位置），點陣的排列有幾種

16、，在 二維模擬中可采用簡單正方點陣或三角點陣，三維模擬大都采用簡單立 方點陣，也有一部分采用面心立方點陣、體心立方點陣或密排六方點 陣。取向?qū)傩钥偣灿蠶個（Q>2），則在模擬時首先對每個原子隨機賦予 一個取向?qū)傩灾礢 （1<S<Q） o Q值一般取32、48、64等值。然后將所有取向取向相同的相鄰原子組合成晶粒，由此構(gòu)造原始的 模擬晶粒組織。Site considered First shell 0 Second shell ： Third shellNN考慮單元的6個最近鄰格點與12個次近鄰格點以及8個第三近鄰的格點oNNE產(chǎn)丿工（1心宀），心宀N NNH工（&冋-

17、1），z=l >1界面能由描述原子相互作用的哈密爾頓算子來定義。下式中J>0可以理解 為相鄰原子間的相互作用能。對于任意格點，其界面能E為：O：Sj 豐 S 1, S： = S jOQ.HSj其中，丿是正的常數(shù)，正比于晶界能，為相鄰單元對晶界能所做的貢獻。 Si、S丿分別對應(yīng)于單元/和/的取向（1WS/WQ）j是Kronecker delta函數(shù)，當&和Q取向相同時，等于1,反之則 等于0。材料設(shè)計優(yōu)化與生物醫(yī)用材料研究室材料設(shè)計優(yōu)化與生物醫(yī)用材料研究室(b)波茨模型中晶界表征示意圖A E=2-4 = -2<0E°=4Ej=25555 55 55 5材料設(shè)計

18、優(yōu)化與生物醫(yī)用材料研究室：«隨機地選取一個單元i，計算該單元的界而能E 0瞬普.個可能的取向中隨機地選取-個賦給該單元,并且計算其自材料設(shè)計優(yōu)化與生物醫(yī)用材料研究室比較單元取向改變前后的能量差A(yù) E(A E = E1-EO)如果A E WO,則新取向被接受，即晶界遷移的驅(qū)動力為晶界能的減小; 否則，新取向以概率W=exp(- A E/kT)被接受。用公式表示，單元成功再取向的概率為：? AE<0 ? AE>0其中，k為Boltzman常數(shù);T為溫度。懿界般騁課髀蠶從阱為M晶界相應(yīng)的向源子原來的疋蒙翥聲堯霽歸整矍窗瘠成鶉黠健黔嚟般航橙賀彖庫界材料設(shè)計優(yōu)化與生物醫(yī)用材料研究室

19、材料設(shè)計優(yōu)化與生物醫(yī)用材料研究室POTTS MODELFRONT TRACKING MODEL材料設(shè)計優(yōu)化與生物醫(yī)用材料研究室2、改進的Monte Carlo算法標準的Monte Cewlo算法在實際應(yīng)用中存在一些不足之處:一是耗時且效率低，尤其在晶粒長大模擬后期，有效取向概率的大大 降低,使晶粒長大指數(shù)遠低于理想值0.5。二是取向總數(shù)Q不易確定。Q值取得過大時，很多再取向嘗試導(dǎo)致系統(tǒng) 能量升高而被拒絕，使晶粒長大過程非常緩慢；如果Q太小，本來 分離的不同晶粒容易“融合"成一個晶粒，出現(xiàn)所謂的“晶粒粗化”現(xiàn) 象，這不是晶粒長大過程的真實現(xiàn)象，小的Q值對于模擬晶界能各 向異性組織和織構(gòu)

20、組織的晶粒長大也會帶來問題。 Radhakrishnan和 Zacharial995年提岀了一種 改進的 Monte Carlo算法：由于只有晶界附近的原子跳到鄰近晶粒中才能產(chǎn)生 晶界的遷移，因此，在單元再定向時,不再考慮所有可能的 取向值Q-1,而僅僅計算該單元周圍鄰近單元的取向值， 這顯然能夠加速晶粒長大的模擬過程，而且可以允許Q 值擴大到無窮大而不影響效率。這種改進算法的時間效率比標準的Monte Carlo算法 提高了大約10倍,晶粒生長指數(shù)接近理論值。改進的方法 已經(jīng)被廣泛接受。(Radliakiisluian B, Zacharial T, Metall. Mater.Trans,

21、 1995?26A(1): 167-180)3、宋曉艷博士改進的Monte Carlo算法對Monte Carlo算法的另一個改進來自于我們課題組， 進一步改進了單元的選取方法，提岀“逐一選取”的方法， 所有單元全部提取并完成一次再取向嘗試記為一個M C So宋曉艷博士指出，對單元的隨機選取徒增時間，既 然系統(tǒng)內(nèi)所有單元均需進行再取向嘗試，這種隨機選取 沒有什么意義，反而增加了程序的運算量，浪費了模擬 時間。這種改進算法使模擬效率又有了很大的提高。采用該 算法，得到了較好的模擬結(jié)果(n=0.49土0.1) oSong Xiaoyan, Liu Guoquan,et al. A modified

22、 Monte Carlo method for grain growth simulation. Progress in Natural Science, 1998, 8(1): 92-97.4、秦湘閣博士改進的Monte Carlo算法本課題組秦湘閣博士吸收了元胞自動機法的思想，進一步 改進了Monte Carlo算法。使在一個MCS內(nèi)，各個單元同 時進行再取向嘗試，每個單元再取向的結(jié)果只取決于此刻 它的臨域狀態(tài)。取向改變的單元緩保存在一個臨時數(shù)組， 當所有單元都進行完再取向嘗試完成之后（即1MCS）, 用臨時數(shù)組中的數(shù)據(jù)更新模擬點陣數(shù)組。秦湘閣不同算法和參數(shù)的晶粒長大P0血模型研究博士學(xué)位

23、論文北京:北 京科技大學(xué)，1997材料設(shè)計優(yōu)化與生物醫(yī)用材料研究室利用動力學(xué)波茨蒙持卡羅覆型對三堆骷校生長的模姒結(jié)果(Mehnert, 1998)應(yīng)用動力學(xué)波茨蒙特卡羅棋型當存在粒子時對二維晶稅生長的橈扭(Doh訥y等人.1990材料設(shè)計優(yōu)化與生物醫(yī)用材料研究室二維復(fù)相材料仿真r =50MCSr=100MCS/=5t)0MCS二維復(fù)相材料顯微組織及其演變過程其中第二相粒子面積分數(shù)為仕0.01 （宋曉艷博士）三維復(fù)相材料仿真仿真的三維復(fù)相材料截面組織與真實組織的比較(a)仿真織織(b)含有穩(wěn)定碳化物粒子的高速鋼的截而組織(劉國權(quán),宋曉艷.自然科學(xué)進展? 1999, 2: 161)材料設(shè)計優(yōu)化與生

24、物醫(yī)用材料研究室5. Monte Carlo模擬得到的一些結(jié)果1. 晶粒組織3D可視化與多角度觀測2. 晶粒長大動力學(xué)3. 晶粒尺寸分布4. 晶粒拓撲特征材料設(shè)計優(yōu)化與生物醫(yī)用材料研究室1啟膘iL番 心占L 盤二穢C 罟二嚨n _1 pnm-mo© 一 miu-sb-jtuy 一uesfto 墨呱肛a 一喲iI辭理DBI述口si 卻 口嗣(創(chuàng)李禪+<* 黛m2ami 9HSEOdbIK鞍雖電無嶇CKS二Ddo系列截面觀測(MCS=3000)作者：王浩個體晶粒的3D可視化與多角度視圖(F = 31)個體晶粒的3D可視化與多角度視圖(F=14)個體晶粒的3D可視化與多角度視圖（F=

25、12）材料設(shè)計優(yōu)化與生物醫(yī)用材料研究室材料設(shè)計優(yōu)化與生物醫(yī)用材料研究室材料設(shè)計優(yōu)化與生物醫(yī)用材料研究室3D晶棱可視化材料設(shè)計優(yōu)化與生物醫(yī)用材料研究室正常晶粒長大動力學(xué)P.A貝克(Beck)等人根據(jù)實驗首次提出在晶粒長大過程中平均晶粒尺寸 與等溫退火時向成拋物線關(guān)系< Rt= Ct"上式一般稱為貝克動力學(xué)關(guān)系或貝克公式，其中c是一個系統(tǒng)常數(shù)， 為動力學(xué)時間指數(shù)。 Burke和Turnbull推導(dǎo)出晶粒長大動力學(xué)的拋物線關(guān)系<Rt>m -<Rq >m=Kt<Rtm = Kt< Rt >= Ktnn般稱為晶粒長大指數(shù)，理論值n=0.5 (不考

26、慮雜質(zhì)原子阻礙作用及 晶界能和晶界遷移率為常數(shù)的情況下)。正常晶粒長大動力學(xué)材料設(shè)計優(yōu)化與生物醫(yī)用材料研究室材料設(shè)計優(yōu)化與生物醫(yī)用材料研究室0500100015002000250030003500Time t(MCS)Simulation dataFitting curve n=0.501648 64vR>=(K(tt°)+vR >')材料設(shè)計優(yōu)化與生物醫(yī)用材料研究室人們對晶粒長大過程中平均晶粒尺寸 人 與時間T 的拋物線動力學(xué)關(guān)系認識已非常深入，然而對單個三 維晶粒的生長速率的研究仍相當不完善。這跟晶粒組 織本身的復(fù)雜性密切相關(guān)一晶粒是由晶界、三晶棱和 四晶交點

27、這些單元組成的復(fù)雜結(jié)合體，晶粒長大需要 這些單元的共同作用才能完成。Russ J C, Delioff R T. Practical Stereology (2nd Edition). Plenum Press, New York, 1999.12.材料設(shè)計優(yōu)化與生物醫(yī)用材料研究室b)c)<DImage by Professor Qin Xiangge晶粒尺寸不僅是反映晶粒長大過程中晶粒組織特征的重要的定量參 數(shù)，而且由于直接影響材料的力學(xué)性能及物理性能，也是實際生產(chǎn) 工藝控制中一個重要的材料微觀組織特征量。粗略定性分析：細晶、粗晶、尺寸大致均勻等半定量分析：各種金相標準(如ASTM標準)定義的晶粒號更精確的定量分析：平均晶粒尺寸R、單位體積中的晶粒數(shù)Nv、 單位面積截面上的晶粒數(shù)Na、單位長度測量 線上的晶粒數(shù)Nl等定量指標。晶粒尺寸分布離散程度：標準方差或者方差系數(shù)。要全面地反映一個試樣中的晶粒尺寸分布情況就需要用尺寸分布函 數(shù) (grain size di