數值模擬計算方法

❶ 數值模擬主要過程和步驟

1、首先要建立反映問題（工程問題、物理問題等）本質的數學模型。

具體說就是要建立反映問題各量之間的微分方程及相應的定解條件。這是數值模擬的出發點。沒有正確完善的數學模型，數值模擬就無從談起。牛頓型流體流動的數學模型就是著名的納維—斯托克斯方程（簡稱方程）及其相應的定解條件。

2、尋求高效率、高准確度的計算方法

由於人們的努力，目前已發展了許多數值計算方法。計算方法不僅包括微分方程的離散化方法及求解方法，還包括貼體坐標的建立，邊界條件的處理等。這些過去被人們忽略或迴避的問題，現在受到越來越多的重視和研究。

3、開始編製程序和進行計算

實踐表明這一部分工作是整個工作的主體，占絕大部分時間。由於求解的問題比較復雜，比如方程就是一個非線性的十分復雜的方程，它的數值求解方法在理論上不夠完善，所以需要通過實驗來加以驗證。正是在這個意義上講，數值模擬又叫數值試驗。應該指出這部分工作決不是輕而易舉的。

(1)數值模擬計算方法擴展閱讀：

數值模擬的發展史：

1955年Peaceman與Rachford研發的交替隱式解法(ADI)是數值模擬技術的重大突破。該解法非常穩定,而且速度快，所以迅速在包括石油,核物理，熱傳導等領域得到廣泛應用。1958年Douglas,Jim和Blair,P.M第一次進行了考慮毛管壓力效果的水驅模擬。

60年代數值模擬技術的發展主要在數值解法，第一個有效的數值模擬解法器是1968年Stone推出的SIP(Strong Implicit Procere)。該解法可以很好地用來模擬非均質油藏和形狀不規則油藏。

Stone在70年代發表了三相相對滲透率模型，由油水和油氣兩相相對滲透率計算油、氣、水三相流動時的相對滲透率，該技術現在還廣為應用。70年代另一項主要成就是Peaceman提出的從網格壓力來確定井底流壓的校正方法。

參考資料來源：網路—數值模擬

❷ 試說明數值模擬方法的特點，它與理論研究，實驗研究有什麼關系

一、數值模擬的特點：

1、藉助於計算機實現。在計算機上實現一個特定的計算，非常類似於履行一個物理實驗。這時分歲兄則析人員已跳出了數學方程的圈子來對待物理現象的發生，就像做一次物理實驗。

2、運算高效且精準。因為數值模擬利用的是計算機系統中穩定的計算。

二、數值模擬與理論研究、實驗研究的關系：

1、數值模擬是塵塵以理論研究為基礎進行的模擬研究，沒有科學的理論支持，無法做到高效的數值模擬。

2、數值模擬類似於實驗研究，是利用計算機的實驗研究。即利用計算機來進行實驗。

(2)數值模擬計算方法擴展閱讀

數值模擬的步驟為：

1、首先要建立反映問題（工程問題、物理問題等）本質的數學模型。具體說就是要建立反映問題各量之間的微分方程及相應的定解條件。這是數值模擬的出發點。

2、數乎棚學模型建立之後，需要解決的問題是尋求高效率、高准確度的計算方法。計算方法不僅包括微分方程的離散化方法及求解方法，還包括貼體坐標的建立，邊界條件的處理等。

3、在確定了計算方法和坐標系後，可以開始編製程序和進行計算。實踐表明這一部分工作是整個工作的主體，占絕大部分時間。

4、在計算工作完成後，大量數據只能通過圖像形象地顯示出來。因此數值的圖像顯示也是一項十分重要的工作。利用錄像機或電影放映機可以顯示動態過程，模擬的水平越來越高，越來越逼真。

❸ 數值方程與數值模擬

常用的數值計算方法有有限差分法和有限單元法。由於有限單元法中的集中儲量有限元方法較通常的有限元法具有更多的優點，而且在邊界條件的處理上，集中儲量有限元法比有限差分法更符合實際，它考慮了邊界節點的均衡單元的儲水量變化（吳金全，1989）。

圖1.4.3 均衡區域示意圖

（一）集中貯量有限元公式推導

取單位水平面積、高度為計算層厚度的土柱進行研究（圖1.4.3），將土柱（計算區域）垂直向上剖分為n個單元，空間步長為Δz，節點編號為0，1，2，…，n-1，n，在Δt時段內（Δt=t^j+1>-t^j），對任一內節點i所代表的均衡區z_i-1/2到z_i+1/2（圖1.4.3）之間的土體列水量均衡方程（暫先不考慮極系吸水項）。

1.內節點（i=1，2，…，n-1）

由達西定律

（h 方程）得：

（1）通過z_i-1/2斷面的水流通量（流入量）為：

土壤水鹽運移數值模擬

（2）通過z_i+1/2斷面的水流通量（流出量）為：

土壤水鹽運移數值模擬

（3）均衡區域Δz內儲水量的變化量（增量）為：

土壤水鹽運移數值模擬

根據質量守恆原理（流入量-流出量=儲存量的變化量）得：

Δq_i=q_i-1/2-q_i+1/2 （1.4.19）

將式（1.4.16）、式（1.4.17）、式（1.4.18）、式（1.4.19）代入得：

土壤水鹽運移數值模擬

式（1.4.20）中，負壓h及參數C和K在時間上取時段末j+1時刻的值，並整理得：

土壤水鹽運移數值模擬

與有限差分方程比較，集中儲量有限元推導出的有限元方程式（1.4.21）與隱式差分方程（h方程）是完全一致的。因此，具有無條件穩定和收斂的優良特性，故選用隱式差分格式對數學模型進行數值離散。若在時間上取時段中間j+1/2時刻的負壓h及參數C、K，則可得出與Crank-Nicolsen差分格式完全一致的方程。

若考慮源匯項根系吸水項S，則式（1.4.22）變為：

土壤水鹽運移數值模擬

令：

，

，r₃=Δt得：

土壤水鹽運移數值模擬

式中：i=1，2，…，n-1。

令

土壤水鹽運移數值模擬

將式（1.4.24）代入式（1.4.23）得：

土壤水鹽運移數值模擬

2.邊界節點的處理

（1）上邊界節點i=0處的方程為：

土壤水鹽運移數值模擬

式中：

為

，R^j+1/2，為時段平均入滲強度，E^j+1/2為時段平均表土蒸發強度。

令（1.4.26）式中：

土壤水鹽運移數值模擬

則（1.4.27）式變為：

土壤水鹽運移數值模擬

（2）下邊界節點i=n為第一類邊界節點，h_n已知，故不需列方程計算，這樣第n-1個方程可簡化為：

土壤水鹽運移數值模擬

式中：

土壤水鹽運移數值模擬

3.方程組

綜合內節點和邊界節點方程，從而得如下代數方程組：

土壤水鹽運移數值模擬

方程組式（1.4.31）中：b₀、c₀、f₀按式（1.4.27）式計算，f_n-1按式（1.4.30）式計算，其餘α_i、b_i、c_i按式（1.4.24）計算。

方程組式（1.4.31）用矩陣表示可簡化為：

［A］［H］^j+1=［F］ （1.4.32）

式中：［ A］為系數矩陣；［ F］為常數項列陣；［ H］^j+1為求解未知量的列陣。

這樣，通過數值方法將描述土壤水分運動的偏微分方程轉化為求解代數方程組的問題。方程組式（1.4.31）系數矩陣元素滿足α_ij=0（當｜i-j｜＞1 時），為三對角方程組，所以，可用「追趕法」求解。

（二）方程的線性化與土壤水分運動參數的取值

系數矩陣［A］中的各元素由時段末（j+1）時刻的土壤水分運動參數給出，常數項列陣［F］中的元素除含有已知時段初j時刻的負壓h外，還含有時段末（j+1）時刻的土壤水分運動參數。然而土壤水分運動參數本身又是負壓h的函數，因而求解方程組原則上說是非線性的。在利用數值方法求解土壤水分運動方程時，必須將方程線性化，使求解方程組成為線性代數方程組。

因迭代法計算的誤差可以控制，求得的結果較逼近實際，而且一般可允許選用較大的時間步長（雷志棟等，1988），故選用迭代法進行線性化。

首先取時段初的參數如

作為時段末參數

的預報值，然後解方程組［ A］［H］^j+1=［F］，求得時段末各節點負壓h的第一次迭代值

，根據

及K-h曲線可求得土壤水分運動參數

的校正值。以此參數的校正值作為下一次計算的預報值，然後解方程組可得時段末各節點負壓h的第二次迭代值

。重復上述步驟，直到前後兩次迭代值，第p-1次和第p次迭代值滿足下式為止：

土壤水鹽運移數值模擬

式中：e迭代誤差為任意給定得正的小數，一般取e=0.01。

參數的取值，一般的說，用三點式或幾何平均的方法效果較好（雷志棟等，1988），計算也不復雜，這里選用幾何平均的方法：

土壤水鹽運移數值模擬

同理，根據達西定律（

）（θ方程）可以推導出集中儲量有限元公式的θ方程，該方程與隱式差分方程θ方程完全一致，若在時間上取時段中間j+1/2時刻的含水率θ及參數D、K，則可得到與Crank-Nicolsen差分格式（中心差分）完全一致的方程。θ方程與h方程的求解方法完全一樣，由於本文所研究的是雙層結構問題，而在分界面處θ不連續、h連續，所以選用h方程進行計算，因此對於θ方程這里就不做推導了。

（三）數值模擬

1.模型驗證

進行數值模擬，首先進行模型驗證。模型驗證時，上邊界條件表達式中的θ₁₀由實際觀測資料給出。根據有作物生長條件下土壤水分運動的基本方程和差分方程，在已知初始條件和邊界條件時，模型驗證可以通過以下步驟進行：①根據實測初始負壓剖面的分布，用三次樣條插值給出各節點上的初始值；②計算蒸發量E；③計算根系吸水層厚度L_r及吸水率S；④根據差分方程計算時段末負壓值h。

模型驗證時，以計算起始時刻的實測負壓剖面（或含水率剖面）作為初始剖面，空間步長選用1cm，根據最底部負壓計測點和中子儀測點，剖面深度為120cm（含水率剖面為130cm），時間步長選用1h，迭代相對誤差e≤0.01。計算中輸入的大量信息，如各節點的初始負壓值、降雨量、降雨日期、水面蒸發強度、根層土壤含水率等均以數據文件的形式提供。由於大田蓋200kg/畝和蓋600kg/畝只進行了含水率觀測，計算時先將含水率θ轉化為負壓h，計算結束後再將負壓h轉化為θ。大田蓋400kg/畝有負壓資料，可直接用負壓h計算。上邊界條件由E/E₀-θ關系給出。數值模擬按覆蓋量（200、400、600kg/畝），分生育階段（400kg/畝）進行，模擬計算在PC機上完成。主要模擬的試驗處理有；擬合曲線見圖1.4.4。

（1）I-2蓋200kg/畝，模擬時間為7月31日至8月30日，共31天。

（2）I-3蓋400kg/畝，模擬時間為苗期：6月25日至7月30日，共26天；拔節：7月21日至8月10，共21天；灌漿-成熟：8月11日至9月17日，共38天。

（3）I-4蓋600kg/畝，模擬時間為7月31日至8月30日，共31天。

2.模型驗證結果及討論

根據描述土壤水分運動的定解問題，通過數值模擬可以得到土壤水分運動的動態過程，並用實測結果對模型進行驗證。如果數學模型能夠描述實際的物理過程，排除隨機因素外，模擬得到的土壤水分動態過程（模擬值）與實際觀測得到的土壤水分動態過程（實測值）應該完全吻合。比較圖1.4.4，從圖中可以看出，模擬值與實測值吻合較好，表明本文提出的考慮秸稈覆蓋有作物生長條件下的模型是可靠的，以上討論的數值方法是可行的。不同覆蓋量、不同生育階段，可以用不同的E/E₀-θ經驗公式來反映覆蓋對水分運移的影響。因此，本文提出的模型和數值方法可以用來模擬秸稈覆蓋條件下田間土壤水分的運動，可對田間土壤水分動態作中短期預報。

圖1.4.4 實測值與模擬值對比圖

3.模型的應用——預報

數值模擬的目的之一就是進行預報，根據氣象部門提供的降雨量及水面蒸發強度等氣象資料，使用驗證過的模型進行田間土壤水分動態預報。本文使用實際發生過的降雨量及水面蒸發強度系列資料進行預報，用實測負壓資料檢驗預報結果。程序的運行見圖1.4.5。檢驗的實測資料選用大田覆蓋400kg/畝的資料，分生育階段（苗期、拔節、灌漿-成熟）進行。從圖1.4.6可以看出，預報值和實測值吻合較好。

圖1.4.5 雙層結構有根系吸水項垂向一維土壤水數值模擬框圖

圖1.4.6 預報值與實測值對比圖（大田蓋400kg/畝）

❹ 數值模擬方法

如前所述，基於鏡質組化學反應動力學的煤熱演化史數值模擬方法經歷了簡單函數關系模式、受熱時間－經驗法模式、反應活化能－溫度函數模式和平行反應化學動力學模式四個發展階段（秦勇等，1995）。其中，目前在化石能源地質界最為常用的是TTI法、LOM法和Tissot法，近年建立起來的EASY%R_o方法也已成為數值模擬的重要發展方向。

（一）時間－溫度指數（TTI）法

這種方法的基本思想由Lopatin（1971）首先提出，後經Waples（1980）根據31個盆地402個樣品的實測資料加以完善，建立起了其數值模擬的方法體系，故又稱為Lopatin-Waples法。TTI法遵循化學反應動力學的基本法則來衡量溫度和時間兩個因素對煤化作用的貢獻，即煤化作用速度隨受熱溫度的增高而增大，溫度每增加10℃煤化反應速度增加一倍。

由此，可將TTI數學模式用連續函數加以表示：

山西南部煤化作用及其古地熱系統：兼論煤化作用的控氣地質機理

T=T_s＋G·D

式中：T——煤層受熱的古地溫溫度，℃；

T_s—－古地表溫度，℃；

G——古地溫梯度，℃/m;

D—煤層埋藏深度，m。

上述連續函數模式不便於常規運算，故在實際工作中通常採用分段積分方式，根據如下方程進行計算：

TTI＝∑γⁿ·△t

式中：γ——溫度效應因子，考慮到溫度每增高10℃煤級量值增加一倍的關系，取r=2;

n——溫度指數沖仿，其大小取決於溫度（表5-1）；

△t——煤層在某一溫度段的受熱時間，Ma。

表5-1TTI數學模式中溫度指數（n）與煤層受熱溫度段的對應關系

註：溫度每增（減）10℃，溫度指數增（減）1。

TTI值本身可以作為煤級或有機質成熟度指標，Waples（1980）也提供了TTI值與鏡質組最大油浸反射率的關系（表5-2）。此外，Wood（1988）通過回歸擬合，將二者之間關系採用下式表示：

lg%R_o＝—0.006（lgTTI）³+0.042（lgTTI）²+0.162（lgTTI）—0.397

表5-2時間溫度指數（TTI）與鏡質組反射率（R_o）的相關關系

應予指出，TTI法盡管幾十年來被化石能源地質界所廣泛採用，但仍存在某些明顯的不足有待於改進，國內外學者為此做過一定努力。例如，溫度每碧猜增加10℃煤化速率加快一倍的假設有其局限性，只能在活化能10～25kcal/mol（相當於20～160℃的溫度范圍，1kcal=4186.8J，下同）區間是可靠的（Magoon,1983）。再如，沉積有機質熱降解的反應活化能是受熱溫度的函數，隨成熟度的增加而逐步加大，不同類型有機質達到相同成熟度所需的活化能也有差異，而TTI模式將整個煤化過程中的反應活化能作為常量，顯然與實際情況有出入。

（二）悔判型有機成熟度水平（LOM）法

該方法由Hood（1975）建立，後經Bostick等（1979）完善而成為一種國內外廣泛採用的數值模擬方法，故又被稱為Hood-Bostick方法。

Hood等考慮到總體反應活化能隨受熱溫度而增高的變化趨勢（18～33kcal/mo1），採用鏡質組油浸反射率標定有機質成熟度水平（LOM）（圖5-1）。在模式中採用了有效受熱時間的概念，即溫度不低於最高受熱溫度15℃范圍內的受熱時間，建立起有機質成熟度、溫度及有效受熱時間之間的相互關系（圖5-2）。

圖5-1有機質成熟度水平（LOM）與其它有機質成熟度指標的關系（引自周中毅，1990）

Ⅰ—褐煤；Ⅱ₁—亞煙煤；Ⅱ₂—高揮發性煙煤；Ⅱ₃—中揮發性煙煤；Ⅱ₄—低揮發性煙煤；Ⅲ₁—半無煙煤；Ⅲ₂—無煙煤

1Btu=1055.06J

在實際工作中，LOM法採用圖解方法求取有機質成熟度或鏡質組反射率（圖5-2）。國內外應用效果表明，利用這種方法推算出的盆地古地溫溫度是比較可靠的，與由綜合研究所得出的結果以及盆地實際情況較為吻合（Vote,1981；周中毅等，1983,1984,1985）。

❺ 數值模擬流程

不同的軟體進行數值模擬時所需的參數、計算方法、剖分格式等不盡相同，數值模擬的過程也不同，但大致相同，本文以TOUGHREACT為例介紹CO₂地質儲存數值模擬的流程。

（一）研究范圍的確定

一般情況下，獨立的天然水文地質系統是計算區最好的選擇，它具有自然邊界，便於較准確地利用其真實的邊界條件，避免人為邊界在資料提供上的困難和誤差。但是在實際工作中，常常不能完全利用自然邊界，這就需要充分利用勘察和長期觀測資料等建立人為邊界。在確定計算區域時，除了保證范圍足夠大以外，還應使假定的邊界條件盡可能接近真實狀態。

計算范圍的劃定應充分考慮研究目的、區域地質構造、儲層岩性、儲層岩石礦物組成及地下水化學成分等多方面因素。數值模擬時間根據研究目的不同稿謹具有不同的時間尺度。就CO₂地質儲存數值模擬而言，如果不考慮地球化學作用，封存系統在1000年數量級的模擬時間內基本上已達到平衡或穩定。在劃定邊界時還應考慮CO₂在儲層中的擴散距離，與研究區地質模型的孔隙度，滲透率等參數關系密切。為了保證所選模型範圍邊界在模擬期內不影響模擬結果，盡量通過具有相同地質條件的天然CO₂氣田（藏）進行類比，確定大體的計算范圍的邊界。如果考慮地球化學反應，由於CO₂注入引發的水－岩－氣反應對圍岩岩性改變較顯著，制約著CO₂注人的速度和徑向運移的距離等。

（二）明確研究目的

在進行數值模擬以前首先要明確利用數值模擬技術要解決什麼樣的問題。對於CO₂地質儲存工程而言，進行數值模擬的目的主要是在CO₂地質儲存工程實施前，通過數值模擬技術對工程的選址、方案設計進行優化，工程實施期技術指導、運沖森行期監測及後期CO₂泄漏的風險評估等進行預測，以指導項目科學、合理地實施，將CO₂泄漏風險降至最低。

研究目的決定著前期資料的收集類型、地質建模的側重點、地質模型離散的精密度以及初始、邊界條件的處理方式等過程。

（三）資料的收集整理

1）通過遙感、綜合地質調查、物探、鑽探和各類樣品測試分析等手段獲取場地深部地層岩性、地質構造、水文地質、水文地球化學、岩石礦物資料和數據；

2）搜集和分析CO₂地質儲存場地地質岩性、區域構造格架、活動斷層與地震活動情況等；

3）採用鑽井岩心、測井和地震反射方法，調查CO₂地質儲存場地目標儲層和蓋層的空間分布形態，埋深、厚度和規模等；

4）使用X射線衍射、掃描電鏡等方法研究分析封存場地岩石礦物組成、孔隙結構特徵及其物理化學性質；

5）通過採取淺部、深部含水層水樣進行水質全分析，獲得儲蓋層地層水及淺部含水層初始水化學成分。

不同的數值模擬軟體其數學模型的數值解法不同，空間離散方式也不盡相同，所需的模型參數也有一定的差異，表9-1即為TOUGHREACT數值模擬所需要的主要參數。

表9-1 CO₂地質儲存模擬過程中需要的主要參數（以TOUGHREACT為例）

圖9-3 網格剖分

網格剖分對計算的精度及計算的效率有很重要的影響。精度越高對模擬結果刻畫的越精細，但是數據的計算量越大，對計算機的要求也越高。建議在進行地質模型剖分時先採用較粗的網格剖分，如果這種剖分方式下模擬結果合理然後再進行精細化剖分，用於對模擬結果更加詳細的刻畫。

2.參數和初始條件

初始條件是指在初始時刻（t=0）時研究區內求解數學模型主要狀態變數的初始值。選擇的應用軟體不同所需的狀態變數數量、種類不同。如TOUGHREACT所需的初始主要狀態變數包括壓力、溫度和組分濃度的空間分布。地質參數包括孔隙度、滲透率、密度、壓力、溫度、毛細壓力等參數值。這些數值一部分採用室內實驗測得，另一部分採用參考文獻的經驗值；地層水的化學成分的初始值採用實際地層水的化學分析，主要是8大離子的濃度、鹽度和pH 等。如果研究區深部地層中的水樣難以獲得，如蓋層，則採用靜態平衡的方法，利用具有與儲層相同鹽度的鹹水與含有原生礦物的地層岩石在原地層環境下進行化學反應，獲取平衡狀態下的地層水化學成分的初始值；通過岩礦分析、電子掃描、Ⅹ衍射等手段，獲得組成CO₂地質儲層蓋層的原生礦物成分體積含量初始值，並根據原生礦物的組成合理判斷次生礦物。

從原則上講，初始時刻是可以任意取定的，只要該時刻所需的參數和狀態變數值已知即可。因此我們不應該把初始條件理解為研究系統的初始狀態。具體如何取，應該視問題的需要、資料來源、計算方便與否等因素而定。

3.邊界條件

邊界條件是某一實際問題數學模型具有定解的必要條件之一。地下水流問題和溶質運移問題邊界條件的定義不盡相同，但一般概化為以下三種。

（1）一類邊界條件（Dirichlet條件）

解決水流問題時，此類邊界條件為在邊界上所有點的水頭是給定的；對於溶質運移問題，一類邊界條件是指研究區邊界上的溶質濃度分布已知。解決CO₂—水兩相流動問題時，此類邊界條件為在邊界上所有點的壓力是給定的。

（2）二類邊界條件（Neumann條件）

當已知某一邊界的單位面積流入或流出的流量時，可視作解決流動問題的二類邊界；相對溶質運移來講，此類邊界又稱給定彌散通量邊界，即邊界上的彌散通量隨時間變化規律已知。

（3）三類邊界（Cauchy條件）

當研究區一部分滿足一類Dirichlet條件，而另一部分滿足二類Neumann條件時，這類問題稱為混合邊界問題，稱為三類邊界。對溶質運移而言，此類邊界為邊界上溶質通量隨時間變化規律已知。

在CO₂地質儲存數值模擬過程中，由於儲層地層多在800m以下，地質模型的頂部和底部根據實際需要可以處理為不透邊界；為了避免邊界對模擬結果的影響，研究區的范圍一般比實際CO₂所能運移到的范圍大得多，因此，在處理四周邊界時一般設置為無窮一類邊界或不透邊界。在確定邊界條件時，應根據水文地質條件以及現有的資料來綜合考慮。

4.源匯項處理

在多孔介質中流動和溶質運移的問題中，對流、水動力彌散和溶質源或／和匯，是決定含水層中任一內點上溶質質量時變率的兩大因素。源匯項問題在水質與水量計算中以及正確處理對流－彌散方程和滲流基本微分方程中佔有重要地位。作為源匯項的方式很多，如越流補給、含水層彈性釋放補給以及抽（注）井的補給等。

對於深部鹹水層CO₂地質儲存系統而言，系統的頂部一般為具有低滲、低孔的泥岩、頁岩等緻密性岩層，越流補給較難發生。整個CO₂地質儲庫系統的源匯項主要指對流（如側向邊界）和抽（注）井。

（八）模型的校正與驗證

模型識別是建立地下流體數值模型最重要的環節之一，正確理解和進行擬合對於提高數值模型的模擬性是至關重要的。在有實測結果的情況下如示範工程，可將模擬結果與實測結果進行比較，對相關參數進行適當合理的調整，使模擬結果在給定的誤差范圍內與實測結果吻合。若誤差較大，應該重新檢驗概念模型的可靠性，甚至重新建立概念模型。在識別校正以後，應採用校正好的模型繼續計算，並與未用來識別校正的實際數據比較，驗證模型的准確性和可靠性。若存在較大誤差，需重復前面的過程。在沒有實測結果的情況下，數值模型的可靠性可通過類比相關資料或根據個人經驗和理論判斷。

（九）模擬預測

模型預測是實施數值模擬技術的主要目的。對於CO₂地質儲存工程而言，由於CO₂地質儲存技術的提出為時尚短，針對CO₂在深部鹹水層中的運移、擴散、與地層水和圍岩產生的化學反應，以及由於CO₂灌注引起的儲蓋層物理、化學性質變化研究均處於研究和發展階段。因此，在工程實施過程中急需具有技術指導性的工具產生，避免造成投資的浪費和CO₂泄漏等風險的出現。

利用經過識別校正與驗證過的數值模型對CO₂地質儲存過程進行模擬預測，有針對性地對模擬數據進行後期處理，如統計分析、比較等手段對結果進行解譯，以此達到場地的優選，目標儲層灌注能力、儲存潛力的評估，CO：擴散運移途徑和速度、不同捕集方式封存量及它們之間的時空轉化等過程的詳細刻畫與模擬模擬等目的。同時可以預測CO₂在已有、重新激活或新生成的裂隙中逃逸的可能性及時間、CO₂泄漏風險評估以及評價CO₂泄漏對淺層地下水的水質、水量及對地表環境的影響等。

上述結果的分析只是數值模擬技術所能解決問題的冰山一角。對於數值模擬結果的處理要根據所研究的目的進行有針對性的提取和解譯。通過對處理後的數據進行總結分析，發現問題從而解決問題，並掌握內在規律，為CO₂地質儲存工程的前期設計、工程實施、中期監測管理提供理論支持和科學的技術指導，並可以提前開展風險預測，盡早制定預案防範CO₂地質儲存工程實施及運行過程中可能出現的隱患。

❻ 湍流的數值模擬方法

現有的湍流數值模擬方法有三種：直接數值模擬、雷洛平均模擬和大渦模擬。

計算復雜程度：雷洛平均模擬-大渦模擬-直接數值模擬

湍流的運動可以看做是時均運動與隨機運動的疊加，按照雷洛時均法，運動中的變數為 dt

不可壓縮時均運動控制方程組之所以出現方程組不封閉，3個速度，物理量，壓力，6個雷洛應力（需求解的位置函數較方程數多），在於方程中出現了湍流脈動值的雷洛應力項，要使得方程組封閉，必須對雷洛應力做出某些假設，即建立應力表達式（或者引入新的湍流方程），通過此表達式把湍流的脈動值與時均值等聯系起來。基於某些假設所得出的湍流控制方程，稱為湍流模型。

基於這些假設有（皮輪世應力：壓應力與切應力）

1、雷洛應力模型：構建雷洛應力補充方程；

2、湍流黏度類模型：引入 湍動黏度 或渦旋黏度，然後湍流應力成為湍動黏度的函數，整個計算在於確定湍動黏度；

3、零方程模型：不使用微分方程，採用代數關系，將湍動黏度與時均值聯系起來；

4、一方程模型：建立了湍動能的運輸方程；

5、兩方程模型：補充兩個微分方程使湍流時均控制方程組封閉的一類處理方法。

湍流 兩方程模型

湍流動能耗散率(turbulent dissipation)是指在分子粘性作用下由湍流動能轉化為分子熱運動動能的速率。通常以單位質量流體在單位時間內損耗的湍流動能來衡量，以ε表示。湍流速度在空間上存在著隨機漲落，從而形成了顯著的速度梯度，在分子粘性力作用下通過內摩擦不斷地將湍流動能轉化為分子運動的動能。大氣湍流的動能耗散主要發生在大小為毫米數量級的湍渦。(%)

湍流動能是湍流速度漲落方差與流體質量乘積的1/2。（J）TKE/M=1/2(ui'^2)

K 越大表明湍流脈動長度和時間尺度越大， ε 越大意味著湍流脈動長度和時間尺度越小，它們是兩個量制約著湍流脈動。

標准 模型（standard model）： 分別建立了 與 的運輸方程， 為湍動能， 為湍動能耗散率，還是用來求解燃肢 。

適應性：（1）模型中的相關系數，主要根據一些特定條件下的試驗結果來確定的；（2）給出的 模型是針對湍流發展非常充分的湍流運動來建立的，即針對 高Re湍流模型 。而當Re較低時（例如近壁面區的流動），湍流發展不充分，湍流的脈動影響可能不如分子黏性影響大，在近壁面可能再現層流。常用的解決壁面流動方法有：一種是壁面函數法；一種是採用低Re的 模型。；（3）標准 用於 強旋流、繞彎曲壁面流動或彎曲流線運動時，會產生一定的失真 。在標准 模型中，對於雷洛應力的各個分量，假定湍流黏度是相同的，即是各向同性的標准，但是在彎曲流線的情況下，湍流時各向異性的。

RNG   模型：通過在大尺度運動項和修正後的黏度項中體現小尺度的影響，而使這些小尺度運動系統地從控制方程中去除。建立的 方程與標准 模型中的方程很相似，與標準的 模型的主要變化有：（1）通過修正的湍動黏度，考慮了平均流動中的旋轉及旋轉流動情況；（2）在 方程中的產生項增加了一項，從而反映了主流時均應變率Eij。這樣，RNG  模型產生項不僅僅與流動情況相關，而且在同一問題中還是空間坐標的函數。 從而可以更好的處理高應變率及流線彎曲程度較大的流動 。RNG  模型仍然是針對充分發展的湍流，而對近壁面區（低Re數）的流動與標准 模型有著相同的困擾。

Realizable  模型：標准模型在應對時均應變率特別大的時候，有可能導致負的正應力，為了使流動符合湍流的物理定律，需要對正應力進行某種數學上的約束，為了保證這種約束的實現，湍流黏度計算式中的系數 應與應變率聯系起來。提出Realizable  模型。與桐備標準的 模型的主要變化有：（1）湍動黏度計算公式發生變化，引入了旋轉和曲率有關的內容；（2） 方程發生大變化；（3） 旋轉均勻剪切流、包含有射流和混合流的自由流動、管道內流動、邊界層流動，以及帶有分離的流動 。

近壁區使用 模型的問題及對策

有固體壁面的充分發展的湍流流動，沿壁面法線的不同距離上，可將流動分為壁面區和核心區兩個部分。核心區為完全湍流區。壁面區分為3個子層：（1）黏性底層；（2）過渡層；（3）對數律層。以黏性力與雷洛切應力的相對大小劃分的。

壁面函數法：將壁面上的物理量與湍流核心區域內待求解的未知量直接聯系。它必須和高Re的 模型配合使用。基本思想為對湍流核心區的流動用 模型求解，而在壁面區不進行求解，直接使用半經驗公式將壁面上的物理量與湍流核心區內求解變數聯系起來。劃分網格時不需要對壁面區加密，只需把第一個內節點布置在對數律成立的區域內，即配置到湍流充分發展的區域。也可以對壁面區網格加密，以得到近壁區物理量分布。

低Re的 模型：為了使得數值計算從低Re區域一直進行到固體壁面上（該處Re=0），主要的方法是在輸送方程中考慮：壁面的黏性、流態的不同、湍流動能的耗散不是各向同性。在使用低Re的 模型進行流動計算時，充分發展的湍流核心區及黏性底層均用同一套公式，這一套公式中當Re數較大時，關於部分修正變成1，且由於黏性底層的速度梯度大因而黏性底層的網格密。 當局部湍流的Re數小於150時，就應該使用低Re的 模型 。

雷洛應力模型

兩方程模型難以考慮旋轉流動及流線曲率變化的影響，直接對Reynolds方程中湍流脈動應力直接建立微分方程，並進行求解，稱之為雷洛應力模型。

RSM由連續性方程、動量守恆方程、應力方程、 方程和 方程構成12個方程組構成封閉的三維湍流流動問題的基本控制方程組，可以通過SIMPLE等演算法求解。如果需要對能量或組分等進行計算，需要建立針對標量型變數（如溫度、組分、濃度）的脈動量的控制方程。

RSM也是針對湍流發展非常充分的湍流運動來建立的，當Re較低時，上述方程不再適用。在近壁區需要採用壁面函數法或低Re的RSM模型。盡管RSM比 模型應用范圍廣，包含更多的物理機理，但是它有很多缺陷。計算實踐表明，RSM雖然能考慮一些各向異性效應，但是不一定比其他模型好，在計算突擴流動分離區和計算湍流運輸各向異性較強的流動時，RSM優於兩方程模型，但是RSM計算量大於兩方程模型。

大渦模擬

大渦模擬的基本假設是：（1）動量、能量、質量及其他標量主要由大渦輸運；（2）流動的幾何和邊界條件決定大渦的特性，而流動特性主要在大渦中體現；（3）小尺度渦旋受幾何和邊界條件影響較小，並且各向同性。大渦模擬過程直接求解大渦，小尺度渦旋用模型來封閉。

LES的控制方程是對N-S方程在波數空間或者物理空間進行過濾得到的。過濾的過程是去掉比過濾寬度或者給定物理量寬度小的渦旋，從而得到大渦漩的控制方程，在變數的方程中加入過濾函數。FLUENT中，大渦模擬只能針對不可壓縮流體的流動。

高雷洛數流動，兩者基本差不多。但是在低雷洛區，重整化群理論的壓網格模型對流動轉捩和近壁流動問題有較好的模擬效果。

轉捩，即從層流到湍流的過渡。流體力學名詞，表徵一種流動現象，英文為transition。轉捩點的計算和預估是設計飛行器的關鍵前提。

❼ 數值模擬的方法有哪些，各自有什麼優缺點，謝謝！

對有條件進行實驗的材料，盡量採用實驗方法，輔以數值模擬檢驗。而在工
程應用中，很多情況下無法進行實驗，如采礦問題等，數值模擬內部程序有相應的計算方法，能模擬較復雜過程。
直觀性與求解速度：實驗直觀性強，數值模擬直觀性不如實驗方法好，較抽象，但可以
快速得到結果。實驗操作復雜。
成本：實驗成本高，數值模擬成本低廉，只需在計算機上進行模擬和數據處理。
施載入荷：數值模擬可以任意施加各種方向的載荷，可以施加實驗方法達不到的條件。
因此數值模擬方法在監測、設備開發、優化、效果預測方面體現了重要價值。
數據採集：實驗只能採集到特定點的的應力應變等數據，不能得到整個材料各點的應力
應變值，而數值模擬方法可以對各個區域、各個測點進行應力分析和位移分析，對實驗進行補充。
數據處理：應將實驗方法和數值模擬方法結合起來使用，分別對結果進行分析後，充分
考慮兩種方法各自的優缺點，互相比較印證，結合理論分析，有針對性地進行數據和結果的修正，才能得到一個比較全面、客觀的結論。
結果可靠性：數值模擬方法在模擬分析過程中，往往要對邊界條件和材料屬性進行簡化，
或多或少對分析結果產生影響，而且結構離散化的形式不同，得到的結果和精度也不同，隨機性比較大，可信度降低。而在實驗中不可避免的客觀、主觀因素也會產生誤差，但是比數值模擬的誤差少得多，可靠性更高。
兩種方法互相檢驗：合理的數值模擬方法對實驗研究和理論分析具有指導作用，可以彌
補實驗工作的不足。實驗與數值模擬結果比較，用來判斷數值模擬方法的可行性。

❽ 有限元數值模擬方法

有限單元法是應用於構造應力場模擬的最廣泛的數學模擬方法。其基本思想是將所研究的地質體以一定的方式（單元形狀和節點個數）簡化為有限個單元組成的離散化模型，再用相應的計算程序求出數值解答。利用有限元法數值模擬，可以利用地質調查和構造解析獲得的較少地質應力狀態的資料來反演區域內各點的應力狀態，從而獲得區域的構造應力場特徵，加深認識區域內的構造演化。目前有限單元法的應用已由彈性力學的平面問題擴展到空間問題、板殼問題，分析對象從彈性材料擴展到塑性、粘彈性、粘塑性和復合材料。

有限元法數值模擬隨著計算機技術的發展在科學計算領域得到廣泛應用，20世紀80年代以來，國際上已有較大型的有限元計算程序達幾百種，其中較著名的有：ANSYS、NASTRAN、ASK、ADINA、SAP等。以ANSYS為代表的數值模擬軟體將有限元分析、計算機圖形學和優化技術相結合，已成為科學計算領域不可缺少的有力工具。

基於本區岩石圈的三維結構特點，我們首先對本區的三層結構相互作用關系進行了模擬。對本區的物理模擬研究，前人已經做過很多工作，其中在對印度板塊擠壓下亞洲中東部的構造模擬中，有的反映出大型走滑斷裂、裂谷和張性盆地以及壓性逆沖斷裂等構造現象，有的反映出多層構造中網路狀流動現象，認為板內變形受塑性流動網路控制（Tapponnier et al.，1982；李建國等，1997）。這些工作往往只反映了本區的某一方面的特性，而無法對本區的構造形態做出動力學的完善解釋。因此在前人的工作基礎上，我們首先建立了本區的一個三層結構模型，其中中上地殼深度根據天然地震資料定為30 km，下地殼以莫霍面為其底界，根據地震測深資料取50 km。因為本模型建立的主要目的是確定岩石圈各圈層之間的作用關系，因此模型底部只考慮到100 km的深度。

❾ 數值模擬

數值模擬（數值法）是對數學模型的一種近似解法，它僅能求出計算域內有限點某個時刻水頭的近似值，這個值在實際應用中可以滿足精度要求。數值法可以解決許多復雜水文地質條件下的滲流計算問題，應用十分廣泛。如用於大中型水源地、地下水的補徑排條件復雜、滲流區形狀不規則、含水介質為非均質各向異性等條件下，確定水頭分布和流量計算。

（一）滲流區域離散化（以二維流為例）

採用數值模擬技術研究地下水的運動，首先將要研究的水文地質模型內的含水層離散化。所謂離散化，就是將要研究的滲流區非均質各向異性含水層，頌滑衡按照一定的方式剖分（分割）成許多相互聯系的小均衡區，在每個小均衡區內是均質各向同性的。在每個小的均衡區內，其含水層參數視為常數；其中心水頭值或有條件下的平均水頭值視為小均衡區內水頭代表值。剖分通常採用兩種形式（矩形、多邊形）進行。

1.矩形均衡域

它是用兩組正交的平行線把均衡區分為許多小的矩形均衡域，如圖7-3所示。在剖分時約定：①定水頭或已知水頭邊界（一類邊界）應從小均衡域的中心通過；②隔水邊界（二類邊界）與小均衡域的邊界重合。這種剖分方法類似於直角坐標系，用適當的編號標定小區域及節（結）點（小均衡域的中心點）。常用的術語有：

圖7-3 滲流區被剖分成矩形小均衡域

（據李俊亭等，1987）

1）點、行、列，點（節點）為小區域的中心點，網格的橫向稱行，豎向稱列。

2）步長，分為空間步長（Δx，Δy，Δz）（圖7-3）和時間步長（Δt）。

3）小區域及節點編號統一記為（i，j），表示小區域及節點位於第i行第j列。

2.多邊形均衡域

由於多邊形均衡域與復雜邊界的幾何形狀比較接近，因此使用較多。它是先按三角形剖分滲流域，再以三角形為基礎構成多邊形均衡域，見圖7-4。常用的術語及注意事項：

1）點元、面元、線元，三角形的邊稱線元，三角形的頂點稱點元（節點或結點），三角形的面積稱面元；

2）要求剖分時三角形的單個內角取30°～90°；

3）滲流區剖後的面積與原面積要吻合，既不要重復也不要開裂。

（二）基本均衡離散方程（以規則網格的有限差分方法為例）

將圖7-3中的（i，j）的均衡區與相鄰均衡域的水量交換關系表示在圖7-5上。

圖7-4 滲流區域三角形

圖7-5（i，j）均衡區的流量關系示意

（據李俊亭等，1987）

1）均衡時段為Δt_n+1：

Δt_n+1=t_n+1-t_n0

表示點（i，j）上t_n時刻的水頭。

2）若（i，j）均衡區內不存在垂向水量交替，則依據水均衡原理有：

地下水動力學

在x軸方向上不同均衡時段分別為：

地下水動力學

式中：

地下水動力學

3）考慮到式（7-14）與式（7-15）的不同，會產生不同的計算結果。計算方案（差分格式）將寫出如下通式：

地下水動力學

式中野做：0≤θ≤1。θ常取3種情況：①當θ=0時稱有限差分法的顯示差分格式；②當θ=1/2時稱有限差分法的對稱（中心）差分格式；③θ=1稱有限差分法的隱式差分格式。

有限差分方程實際上是基本微分方程的近似表達式，其近似程度可用泰勒級數進行分析。通過微分方程的差分表達式，可以看出在利用差分格式代替微分式時，是存在誤差的，即用有限差分方程組模擬地下讓薯水流系統會產生誤差。

（三）對於邊界條件和垂向水量交換的處理

不論是已知水頭的一類邊界或已知流量的二類邊界，計算點落在邊界上，該點就不需要列入均衡方程。垂向水量交換的處理也是如此，若點與抽水井重合，該點已列入均衡離散方程時，抽水量就直接參與該點所在均衡區的水均衡。

（四）均衡離散方程的解算

顯然，在含水層參數和邊界條件都給定的條件下，只要知道某時刻流場中所有點的水頭值，就可計算出下個時間步長的所有點的水頭。即在已知初始條件的基礎上，可以計算不同時刻各點水頭值、不同時刻的流場。對於這類問題的求解方法，從廣泛使用微機處理的角度來看，超鬆弛迭代為許多研究者所採用。

（五）應用

綜上所述，在已知初始條件、邊界條件、垂向水量交換以及給定含水層參數的情況下，可計算滲流區內不同時刻、不同節點的水頭值。當前，不論是在地下水資源評價的水量計算中，還是在礦山開采地下水的疏干計算或在因大面積地下水位下降引起的地質災害防治中，數值法都得到了廣泛應用。

目前有許多地下水數值法計算軟體，適應性強、有較高的模擬性，廣為採用，例如，MOP-FLOW（孔隙水三維有限差分法數值計算軟體），GWMS-3D（二維或三維地下水流和污染物質運移數值模擬軟體）等。

（六）實例

通過實例的學習，使同學們對用數值法求解過程有所了解。這個過程包括：①水文地質條件概化，建立概念模型；②根據水文地質概念模型，建立數值模型；③剖分計算區，整理計算資料；④校正數值模型；⑤驗證數值模型；⑥運用模型進行預報。

實例位於太行山東麓沖洪積扇的交界處。含水層為第四紀鬆散層，上部為細砂和粉砂層，下部為砂卵礫石、粗砂礫石加土層、含粘土礫石層等。上部含水層地下水已被疏干，當前開采層埋深為40～80m，水位埋深多在10m以下，漏斗中心區已達30m。邊沿部分地區水位埋深為2～10m。

1.水文地質概念模型

①含水層底板為隔水粘土層；②含水層主要為非均質各向同性的潛水含水層；③計算區的邊界三面為已知水頭的一類邊界，另一面為不同程度的弱透水層，計算區面積近600km²；④區內有開采井；⑤地下水流為非穩定平面流，水流符合達西流。

2.數值模型

1）微分方程：

地下水動力學

2）初始條件：H（x，y，0）=H₀（x，y）

3）一類邊界條件：H（x，y，t）｜_Γ1=H₁（x，y，t）

4）二類邊界條件：

地下水動力學

式中：W為匯源項，由降水入滲量和井的開采量代數和求出；n為內法線；其他符號同前。

3.剖分計算區並整理計算資料

將計算區剖分為506個小區、230個節點，其中第一類邊界點40個，二類14個，取旱季為模型校正時段，給出10個分區參數並經過試驗給出參數初值。

4.校正數值模型

校正結果表明，微分方程和邊界條件吻合。

5.驗證數值模型

取雨季水位資料，分7個時段進行水位驗證。根據驗證資料繪制高低水位擬合圖以及其他所需擬合圖件，證明擬合程度良好，符合規范要求。

6.模型使用

利用驗證過的符合實際的模擬模型，根據設計水位預計開采量，或根據設計的開采量預計不同時段的水位降低，尤其是漏斗中心的水位降低。

❿ 應力場數值模擬方法

近30年來，人們採用現場測試、實驗室試驗、理論分析與模型試驗等多種方法，使岩土力學研究取得很大進展^{［162～166］}。如今隨著計算機技術的快速發展，岩土力學的研究進入了一個新的階段，其中數值計算方法已成為解決岩土力學問題的重要手段之一。

6.1.1 概述

許多工程分析問題，如固體力學中的位移場和應力場分布分析、電磁學中的電磁場分析、振動特性分析、傳熱學中的溫度場分析以及流體力學中的流場分布等，都可以通過在給定邊界條件下對其控制方程進行求解得到，但是利用解析方法只能求出一些方程性質比較簡單且幾何邊界相當規則的極少數問題。對於大多數實際工程技術問題，由於物體的幾何形狀比較復雜或者問題的某些特性是非線性的，因而一般無解析解。為了解決此類問題，一般採用兩種處理方法：一種是進行簡化處理，將方程和邊界條件簡化為能夠處理的問題，從而得到在簡化情況下的解，但這種方法應用非常有限，且假設過多將會導致錯誤的解；另一種是在廣泛接收現代數學和力學理論的基礎上，藉助於計算機和計算軟體來獲得工程上要求的數值解，這就是目前應用非常廣泛的數值模擬方法。

目前在工程技術領域內常用的數值分析方法包括：有限單元法、邊界元法、離散單元法以及有限差分法。最初常用的是有限差分法，它可以處理一些相當復雜的問題。但對於幾何形狀復雜的邊界條件，其解的精度受到影響。20世紀60年代出現並得到廣泛應用的有限單元法，使經典力學解析方法難以解決的工程力學問題都可以用有限元方法求解。它將連續的求解域離散為一組有限個單元的組合體，解析地模擬或逼近求解區域。由於單元能按各種不同的聯結方式組合在一起，且單元本身又可有不同的幾何形狀，所以能適應幾何形狀復雜的求解域。但有限單元法需要的存貯容量常非常巨大，甚至大得無法計算。由於相鄰界面上只能位移協調，對於奇異性問題（應力出現間斷）的處理比較麻煩，這是有限單元法的不足之處。70年代末期，出現了另一種重要的數值方法為邊界元法。邊界元方法是把求解區域的邊界剖分為若干個單元，將求解簡化為求單元結點上的函數值，通過求解一組線性代數方程實現求解積分方程。上述兩種數值方法的主要區別在於，邊界元法是「邊界」方法，而有限元法是「區域」方法，它們都是針對連續介質，只能獲得某一荷載或邊界條件下的穩定解。對於具有明顯塑性應變軟化特性和剪切膨脹特性的岩體，無法對其大變形過程中所表現出來的幾何非線性和物理非線性進行模擬，這就使得人們去尋求適合模擬節理岩體運動變形特性的有效數值方法。

1971年Cundall，P.A^［167］提出了一種不連續介質數值分析模型——離散單元法。該方法優點在於適用於模擬節理系統或離散顆粒組合體在准靜態或動態條件下的變形過程。離散單元法的基本原理不同於基於最小總勢能變分原理的有限單元法，也不同於基於Betti互等定理的邊界單元法，而是建立在牛頓第二運動定律基礎上。最初的離散元法是基於剛性體的假設，由於沒有考慮岩塊自身的變形，在模擬高應力狀態或軟弱、破碎岩體時，不能反映岩塊自身變形的特徵，使計算結果與實際情況產生較大出入。Maini，T.，Cundall，P.A.^{［168～169］}等人針對剛體單元沒有考慮岩塊自身變形的缺點，利用差分方法提出了考慮岩石自身變形的改進的離散單元法，編制了通用的離散元程序UDEC（Universal Discrete Element Code），將離散元推廣到模擬岩體破碎和變形情況，推動了離散元的進一步發展。我國學者也相繼開展這方面的研究，王泳嘉教授^［170］等將離散單元法應用於采礦工程方面的研究。

6.1.2 FLAC數值模擬方法

（1）概述

數值模擬技術通過計算機程序在工程中得到廣泛的應用。一直到20世紀80年代初期，國際上較大型的面向工程的通用程序有：ANSYS、NASTRAN、FLAC、UNDEC、ASKS以及ADINA等程序。它們功能越來越完善，不僅包含多種條件下的有限元分析程序，而且帶有功能強大的前、後處理程序。

連續介質快速拉格朗日差分法（Fast Lagrangian Analysis of Continua，簡寫FLAC）是近年來逐步成熟完善起來的一種新型數值分析方法。把拉格朗日法移植到固體力學中，即將所研究的區域劃分為網格，節點相當於流體質點，然後按照時步用拉格朗日方法來研究網格節點的運動，這就是固體力學變形研究中的拉格朗日數值研究方法。

FLAC與基本離散元法相似，但它克服了離散元法的缺陷，吸取了有限元法適用於各種材料模型及邊界條件的非規則區域連續問題解的優點。FLAC所採用的動態鬆弛法求解，不需要形成耗機時量較大的整體剛度矩陣，佔用計算機內存少，利於在微機的工程問題。同時，FLAC還應用了節點位移連續的條件，可以對連續介質進行大變形分析。

（2）數學模型

顯式有限差分法的基本方程主要包括：平衡方程、幾何方程、物理方程和邊界條件。在FLAC^3D2.0中採用的拉格朗日描述方程，一般規定介質中一點由向量分量x_i，u_i，v_i，dv_i/dt（i=1，2，3）來表徵，其分別代表位置、位移、速度和加速度分量。

其基本原理和基本公式簡單敘述如下：

空間導數的有限差分近似

三維FLAC方法中採用了混合離散方法，區域被劃分為常應變六面體單元的集合體；而在計算過程中，又將每個六面體分為常應變四面體，變數均在四面體上進行計算，六面體單元的應力、應變取值為其四面體的體積加權平均。

如圖6.1所示，所研究區域任一四面體，節點編號為1～4，規定與節點n相對的面為第n面，設定其內任一點的速度分量為v_i，則由高斯散度定理得

煤岩動力災害力電耦合

式中：V——四面體體積，m³；S——四面體外表面，m²；n_j——外表面單位法向向量分量。

圖6.1 四面體

對於常應變單元，n_j在每個面上為常量，因此通過上式積分可得

煤岩動力災害力電耦合

式中上標f表示f面的變數值，對於為線性分布的速率分量，速度分量的平均值為

煤岩動力災害力電耦合

式中上標l表示節點l的變數值。將（6.3）式代入（6.2）式可得

煤岩動力災害力電耦合

經過變換可得節點速率計算公式：

煤岩動力災害力電耦合

1）平衡方程（運動方程）

顯式有限差分法採用的平衡方程就是人們熟知的牛頓第二運動定律，即

煤岩動力災害力電耦合

式中：F_i——節點合力在i方向分力，N；m_i——節點質量，kg；a_i——節點加速度在i方向分量，m/s²。

作用於各個節點的合力：外力（集中力、均布力、重力等）和內力（單元變形引起的應力在單元節點上的分量）。節點質量是根據節點相鄰單元的面積（體積）和密度，按照面積（體積）加權求出。

FLAC^3D以節點為計算對象，將力和質量均集中在節點上，然後通過運動方程在時域內進行求解。節點運動方程可以表示為如下形式：

煤岩動力災害力電耦合

式中：（t）———t時刻l節點在i方向的不平衡力分量，可以由虛功原理導出；m^l———l節點的集中質量，在分析靜態問題時，採用虛擬質量；而在分析動態問題時，則採用實際的集中質量。

將（6.7）式左端用中心差分來近似，則可得

煤岩動力災害力電耦合

2）變形協調方程——幾何方程

作為連續介質力學，變形體之間必須滿足變形協調方程（幾何方程），否則變形體就會出現分離或嵌入。變形協調方程反映了位移與應變間的關系，對於某一時步的單元應變增量可由下式確定：

煤岩動力災害力電耦合

求出應變增量後，即可由本構方程得到應力增量，各時步的應力增量疊加即可得到總應力，在大變形時，還需根據本時步單元的轉角對本時步前的總應力進行旋轉修正，然後即可由虛功原理求出下一時步的節點不平衡力，進入下一時步的計算。

3）物理方程——本構關系

物理方程反映應力與應變之間的關系，在程序中通常被稱為材料模式或材料模型。在FLAC^3D2.0中提供了10種基本材料模型，它們是：①Null；②Elastic，isotropic；③Elastic，transversely isotropic；④Druck-Prager plasticity；⑤Mohr-Coulomb plasticity；⑥Ubiquitous joint plasticity；⑦Strain-hardening／softening Mohr-Coulomb plasticity；⑧bilinear strain-hardening/softening ubiquitous-joint plasticity；⑨Modified Cam-clay plasticity 和⑩elastic，orthotropic。

本文進行應力場數值模擬時採用的是Mohr-Coulomb應變硬化軟化破壞准則，在FLAC^3D2.0中，Mohr-Coulomb 模型的破壞准則以主應力σ₁，σ₂，σ₃來描述，相應的應變為三個主應變ε₁，ε₂，ε₃。根據Hooke定律，應力、應變增量具有如下表達形式：

煤岩動力災害力電耦合

式中α₁，α₂為材料常數，可以由體積模量K和剪切模量G確定：

煤岩動力災害力電耦合

不失一般性，令σ₁≥σ₂≥σ₃，摩爾—庫侖准則為

其中：

煤岩動力災害力電耦合

式中C，φ分別為煤岩的粘聚力和內摩擦角。

FLAC^3D2.0的Mohr-Coulomb 破壞准則如圖6.2所示。

圖6.2 FLAC^3D的Mohr-Coulomb 破壞准則

本著作中就是選用上述的Strain-hardening／softening Mohr-Coulomb plasticity模型，對單軸壓縮煤岩以及礦山地下煤岩獨巷掘進時圍岩的變形破壞過程進行模擬。

4）阻尼力

對於靜態問題，FLAC^3D2.0在式（6.7）的不平衡力中加入了非黏性阻尼，以使系統的振動逐漸衰減直至達到平衡狀態（即不平衡力接近零），此時節點運動方程變為：

煤岩動力災害力電耦合

式中阻尼力（t）由下式確定：

煤岩動力災害力電耦合

上式中α為阻尼系數，其默認值為0.8；而：

煤岩動力災害力電耦合

5）初始條件與邊界條件

邊界條件包括面積力、集中載荷等應力邊界條件和位移邊界條件。此外也可載入體力和初始應力。在編寫程序代碼時，一般所有的應力和節點速度初始化為零，然後指定初始化應力。集中載荷則載入在面節點上，位移邊界條件則以運動方程形式施加到相應的邊界節點上。

邊界條件分為應力邊界條件和位移邊界條件，應力邊界條件為：

煤岩動力災害力電耦合

式中：F_i———作用於節點i上的力；——作用於邊界上的應力；n_j———邊界上的法線沿j方向的矢量大小；Δs———邊界的長度。

若是位移邊界條件，應將邊界條件以運動方程的形式施加到相應的邊界節點上。

FLAC^3D2.0^［171］與FLAC^2D3.3也是由美國Itasca Consulting Group Inc開發的三維顯式有限差分法程序，它可以模擬岩土或其他材料的三維力學行為。FLAC^3D2.0的計算循環過程如圖6.3所示。

圖6.3 FLAC^3D2.0的計算循環

6.1.3 FLAC數值模擬方法在采礦工程中的應用^{［172～179］}

采礦過程中圍岩活動規律及巷道圍岩穩定性問題涉及岩體力學特性、圍岩壓力、支護圍岩相互作用關系及巷道與工作面時空關系等一系列復雜力學問題。隨著我國經濟建設的高速發展，岩土工程穩定性分析問題日益突出，除采礦工程外，在水利、交通（鐵道和公路）、高層建築的地基等行業也都存在著大量的岩土力學數值計算分析問題。能否用計算機數值模擬分析采礦岩層控制問題和岩土工程問題已成為一個大學岩層控制技術和岩土力學學科水平高低的標志之一。

與ANSYS、ADINA相比，FLAC 和UDEC的最大特點是計算分析岩土工程中的物理不穩定問題，因而特別適用於岩土工程中幾何和物理高度非線性問題的穩定性分析，如采場的采動影響規律，軟岩巷道的大變形問題，采動後的地表沉陷，露天礦的邊坡穩定，水壩的穩定性等問題。

從力學計算方法上講其主要特點

1）可以直接計算非線性本構關系；

2）物理上的不穩定問題不會引起數值計算的不穩定；

3）開放式程序設計（FISH），用戶可以根據需要自己設計程序；

4）既可以分析連續體問題（FLAC），也可以分析非連續體問題（UDEC）；

5）可以模擬分析很大的工程問題；

6）高度非線性問題不增加計算時間。

在采礦工程中，許多學者利用FLAC軟體對采礦過程中圍岩活動規律及巷道圍岩穩定性問題涉及到岩體力學特性、圍岩壓力、支護圍岩相互作用關系及巷道與工作面的時空關系等一系列復雜的力學問題進行了一系列的研究，取得了顯著的效果。梅松華等以施工期監測結果為基礎，在正交設計原理的基礎上，選定反演參數與水平，採用二維顯式差分法FLAC進行彈塑性位移反分析。朱建明等在分析FLAC有限差分程序的基礎上，提出了變彈性模量方法模擬時間因素對巷道圍岩穩定性影響的衰減曲線，為揭示巷道圍岩變形機理和有效指導圍岩支護提供了有效的分析方法。來興平等探討了岩石力學非線性計算軟體FLAC^2D3.3在地下巷道離層破壞數值計算中的應用。康紅普對回採巷道錨桿支護影響因素進行了FLAC分析，認為FLAC^2D3.3在分析幾何非線性和大變形問題方面性能優越。

在煤岩動力災害預測中，這些方法的優點

1）可以提前知道煤與瓦斯突出、沖擊礦壓等煤岩動力災害防治的重點區域；

2）可以得到大范圍內的空間信息；

3）可以提前預測預報煤岩動力災害的危險性；

4）可以確定在採掘過程中，應力的分布狀況和集中程度。

在煤岩動力災害預測中，這些方法也具有以下缺點

1）對實際問題均進行了簡化處理；

2）對於煤岩體的力學特性，如彈性模量、泊松比等力學參數，也進行了簡化，沒有考慮其局部非均質性和各向異性；

3）只能作為一種近似方法使用。