如何用數值模擬的方法計算近似解

1. 數值模擬計算結果與分析

在三層模型計算中，共劃分節點1659 個，劃分單元2316 個。圖7-5a、b 為三層模型結構和約束狀態圖。其中選取模型下襲知底部為南，頂部為正北方向，圖中的實體用深淺不同顏色來區分不同的材料性質。底部的箭頭方向為施加的壓應力方向。

圖7-15 第三個模型不同載入階段應力分布計算結果

圖7-16 第三個模型位移矢量分布計算結果

由於塑性變形與載入時間有關，從不同階段的 y 方向應力分布圖中（圖7-15）可以看出，在載入過程中應力分布在中部出現低值區，反映東西兩側應力、應變大於中部盆地發育區的特點，中部地塊相對兩側向南擠出。

從模型位移矢量分布圖中也可以明顯看出這一特徵，在模型頂部位移矢量圖中（圖7-16a）可以看到模型整體有向北運動的趨勢，但在 EW向受到擾動，在 SN向斷裂「蜂腰」處最為強烈。在模型側面位移矢量圖中（圖 7-16b）可以看到位移隨深度存在明顯變化，越趨於上部，位移越弱。

2. 應力場數值模擬方法

近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）只能作為一種近似方法使用。

3. 數值模擬主要過程和步驟

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

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

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

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

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

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

(3)如何用數值模擬的方法計算近似解擴展閱讀：

數值模擬的發展史：

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

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

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

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

4. 計算力學的數值模擬

力學現象的數學模擬，常常歸結為求解常微分方程、偏微分方程、積分方程、或代數方程。求解這些方程的方法有兩類：一類是求分析解，即以公式表示的解；另一類是求數值解，即以成批數字表示的解。很多力學問題相當復雜，特別是復雜的偏微分方程組，一般難以得出它們的分析解，而用數值方法求解則運算步驟繁復，耗用人力很多，因此在電子計算機出現以前，非不得已不用。20世紀50年代以來，出現了配有現代程序設計語言的通用數字計算機。計算機的快速運算和大存貯量，使解復雜的力學問題成為可能。三十多年來，隨著計算機的改進，數值方法得到廣泛的應用和很大的發展；主要是考慮算得更快、更准、省錢，並為原先不能算的問題構造演算法。數值方法很多，求解偏微分方程數值解，以有限差分方法和有限元法使用最廣；此外，還有變分方法、直線法、特徵線法和譜方法，等等。這些方法的實質絕大多數是將偏微分方程問題化成代數問題，然後再用計算機求未知函大孝數數的數值解。 有簡單、靈活和通用性強等特點。用差分方法求數值解時，須先將自變數的定義域「離散化」，即只企圖算自變數定義域中有限個點的未知函數的近似值。如果自變數只有一個，則可把要計算的區間離散成個線段。如果自變數有兩個，而計算區域是圖1[二變數區域的離散化]所示的矩形，則最簡單的離散方式是把區域分成乘個小矩形。小矩形的長 和寬分別叫作方向和方向的步長。微分方程中出現的偏導數(，)， 在微積分中是差商的極限，在有限差分方法中則代以差商。如圖1[二變數區域的離散化]中點的有的情形可代以差商(()-())/2，有的情形可代以(()-())/，如果有二階偏導滾首數，常常可代以二階差商(()-2()+())/2，其中()、()和()分別表示相應點的值。 如以適當的差商來代替微分方程每一個導數，就得到對應於
原微分方程的差分方程怎樣選差商至關重要。此外，偏微分方程總還要附加邊界或初始條件，這些條件也要用差分形式表示。這樣，對於每個網格點的未知函數慎罩值作出未知量的代數方程組。如果網格分得較密，即步長和都比較小，或與 的數值都比較大，則所得代數方程組的未知量的數目將很大，但藉助計算機，還是可以很快求出解來。由於步長無法取為零，因此用差分方法只能求得原微分方程的近似解。但只要選擇合理的差商和步長，計算結果仍能令人滿意，有時還能得到精度很高的解。有限元法
這種方法是把計算區域剖分成大小不等的三角形（或其他形狀的）單元，然後在各單元上用適當的插值函數來代替未知函數。根據變分原理，可將偏微分方程化成代數方程來求解。這種方法具有很廣泛的適應性，特別適於求解具有復雜邊界形狀和物理條件的問題，而且很容易在計算機上實現。1970年以來已研究出一些適用於廣泛的線性問題的有限元通用程序，對工程設計起很大作用。按照有限元法剖分的思想，把汽車外殼剖分成大小不等的許多三角形單元，而對彎曲邊界只須裁彎取直即可。在應力變化劇烈和要求精確計算的地方，須把單元取得小些；在變化不劇烈的地方則可取得大些。用這種方法不僅可以適應復雜的區域，還可以盡量減少總的單元數目，從而減少未知量的數目。如果在有限差分方法中用矩形網格，則較難處理如此復雜的區域。

5. 數值模擬流程

不同的軟體進行數值模擬時所需的參數、計算方法、剖分格式等不盡相同，數值模擬的過程也不同，但大致相同，本文以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₂地質儲存工程實施及運行過程中可能出現的隱患。

6. 數值模擬

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

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

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

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.模型使用

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

7. fluent 數值模擬計算問題

意思是在有一些網格內，湍流粘性比被限制在10000內，就是說超過了設定值，一般這種情況，兩種解決思路，一是改善網格，一是調整湍流參數，具體情況具體分析啊。

一般是網格問題，或者邊界條件也要再檢查檢查。

不是的話，這里有人寫了很好的解決思路，你可以參考一下，這一般用來解決很麻煩的模擬

這是一個常見的問題，很多人已經問過很多遍了，下面我就描述一下我解決問題的方法：
首先，出現這個警告的最基本的原因是由於邊界條件的設置錯誤，所以如果你確認問題的設定沒有問題，就接著下褲兄面的分析吧。
問題的根源在於求解器在計算兩方程模型的湍動能k和耗散率e或者omega時發生錯誤，因而在計算湍流粘度的時候也會出錯。理想狀況下，當求解收斂時，警告會逐漸消失。但是一般狀況下，結果不會向我們想像的畝純高那樣樂觀，原因就是我們的算例很大，收斂比較困難，並且隨著湍流參數的計算錯誤，狀況會進一步惡化，那麼我們該如何補救呢？
傳統的辦法就是選擇耦合式求解器並啟用，通常會解決問題。但我個人的觀點是如何你的算例是屬於不可壓縮，耦合式求解器並無法有效地解決問題。但我們還有一個更為穩妥的辦法，那就是選擇FAS，增加庫朗數，使粗化級別為4.幾乎所有問題都會收斂的。
當然，如果非要使用分離式求解器，我們可做如下選擇：
如果我只想得到流場的湍流參數的近似值的話，可以先關閉動量方程，只求解k,e或者k,omega知道得到一個穩定的湍流場，這時再開啟所有方程進行迭代求解。這個辦法很有效，但有一個問題，如果網格很龐大，比如上百萬個，那麼光是得到收斂的湍流場就需要一兩天的時間迅尺，那麼我們又該怎麼辦呢？
這時候，我可以准備一個相應的粗網格，邊界條件都一樣，首先在粗網格上得到收斂解，這可能只需要一兩個小時，然後選擇file->interpolate,將得到的數據寫進相應的區域。之後關閉動量計算，只計算湍流量，計算幾步之後開啟所有的求解器繼續迭代，你將會得到所有的收斂解。

最後，這都是我的個人理解，如有誤，請指教。

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

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

9. 數值模擬方法

如前所述，基於鏡質組化學反應動力學的煤熱演化史數值模擬方法經歷了簡單函數關系模式、受熱時間－經驗法模式、反應活化能－溫度函數模式和平行反應化學動力學模式四個發展階段（秦勇等，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）。

10. 數值模擬的計算機方法

有限差分方法(FDM)是計算機數值模擬最早採用的方法，至今仍被廣泛運用。該方法將 求解域劃分為差分網格，用有限個網格節點代替連續的求解域。有限差分法以Taylor級 數展開等方法，把控制方程中的導數用網格節點上的函數值的差商代替進行離散，從而 建立以網格節點上的值為未知數的代數方程組。該方法是一種直接將微分問題變為代數 問題的近似數值解法，數學概念直觀，表達簡單，是發展較早且比較成熟的數值方法。 對於有限差分格式，從格式的精度來劃分，有一階格式、二階格式和高階格式。從差分 的空間形式來考慮，可分為中心格式和逆風格式。考慮時間因子的影響，差分格式還可 以分為顯格式、隱格式、顯隱交替格式等。目前常見的差分格式，主要是上述幾種形式 的組合，不同的組合構成不同的差分格式。差分方法主要適用於有結構網格，網格的步 長一般根據實際地形的情況和柯朗穩定條件來決定。

構造差分的方法有多種形式，目前主要採用的是泰勒級數展開方法。其基本的差分表達 式主要有三種形式：一階向前差分、一階向後差分、一階中心差分和二階中心差分等， 其中前兩種格式為一階計算精度，後兩種格式為二階計算精度。通過對時間和空間這幾 種不同差分格式的組合，可以組合成不同的差分計算格式。

有限元方法的基礎是變分原理和加權餘量法，其基本求解思想是把計算域劃分為有限個互不重疊的單元，在每個單元內，選擇一些合適的節點作為求解函數的插值點，將微分 方程中的變數改寫成由各變數或其導數的節點值與所選用的插值函數組成的線性表達式 ，藉助於變分原理或加權餘量法，將微分方程離散求解。採用不同的權函數和插值函數形式，便構成不同的有限元方法。有限元方法最早應用於結構力學，後來隨著計算機的發展慢慢用於流體力學的數值模擬。在有限元方法中，把計算域離散剖分為有限個互不重疊且相互連接的單元，在每個單元內選擇基函數，用單元基函數的線形組合來逼近單元中的真解，整個計算域上總體的基函數可以看為由每個單元基函數組成的，則整個計算域內的解可以看作是由所有單元上的近似解構成。在河道數值模擬中，常見的有限元計算方法是由變分法和加權餘量法發展而來的里茲法和伽遼金法、最小二乘法等。根據所採用的權函數和插值函數的不同，有限元方法也分為多種計算格式。從權函數的選擇來說，有配置法、矩量法、最小二乘法和伽遼金法，從計算單元網格的形狀來劃分，有三角形網格、四邊形網格和多邊形 網格，從插值函數的精度來劃分，又分為線性插值函數和高次插值函數等。不同的組合 同樣構成不同的有限元計算格式。對於權函數，伽遼金(Galerkin)法是將權函數取為逼近函數中的基函數 ；最小二乘法是令權函數等於餘量本身，而內積的極小值則為對代求系數的平方誤差最小；在配置法中，先在計算域 內選取N個配置點 。令近似解在選定的N個配置點上嚴格滿足微分方程，即在配置點上令方程餘量為0。插值函數一般由不同次冪的多項式組成，但也有採用三角函數或指數函數組成的乘積表示，但最常用的多項式插值函數。有限元插值函數分為兩大類，一類只要求插值多項式本身在插值點取已知值，稱為拉格朗日(Lagrange)多項式插值；另一種不僅要求插值多項式本身，還要求它的導數值在插值點取已知值，稱為哈密特(Hermite)多項式插值。單元坐標有笛卡爾直角坐標系和無因次自然坐標，有對稱和不對稱等。常採用的無因次坐標是一種局部坐標系，它的定義取決於單元的幾何形狀，一維看作長度比，二維看作面積比，三維看作體積比。在二維有限元中，三角形單元應用的最早，近來四邊形等參元的應用也越來越廣。對於二維三角形和四邊形電源單元，常採用的插值函數為有Lagrange插值直角坐標系中的線性插值函數及二階或更高階插值函數、面積坐標系中的線性插值函數、二階或更高階插值函數等。

對於有限元方法，其基本思路和解題步驟可歸納為
(1)建立積分方程，根據變分原理或方程餘量與權函數正交化原理，建立與微分方程初邊值問題等價的積分表達式，這是有限元法的出發點。
(2)區域單元剖分，根據求解區域的形狀及實際問題的物理特點，將區域剖分為若干相互連接、不重疊的單元。區域單元劃分是採用有限元方法的前期准備工作，這部分工作量比較大，除了給計算單元和節點進行編號和確定相互之間的關系之外，還要表示節點的位置坐標，同時還需要列出自然邊界和本質邊界的節點序號和相應的邊界值。
(3)確定單元基函數，根據單元中節點數目及對近似解精度的要求，選擇滿足一定插值條 件的插值函數作為單元基函數。有限元方法中的基函數是在單元中選取的，由於各單元 具有規則的幾何形狀，在選取基函數時可遵循一定的法則。
(4)單元分析：將各個單元中的求解函數用單元基函數的線性組合表達式進行逼近；再將 近似函數代入積分方程，並對單元區域進行積分，可獲得含有待定系數(即單元中各節點 的參數值)的代數方程組，稱為單元有限元方程。
(5)總體合成：在得出單元有限元方程之後，將區域中所有單元有限元方程按一定法則進 行累加，形成總體有限元方程。
(6)邊界條件的處理：一般邊界條件有三種形式，分為本質邊界條件(狄里克雷邊界條件 )、自然邊界條件(黎曼邊界條件)、混合邊界條件(柯西邊界條件)。對於自然邊界條件， 一般在積分表達式中可自動得到滿足。對於本質邊界條件和混合邊界條件，需按一定法 則對總體有限元方程進行修正滿足。
(7)解有限元方程：根據邊界條件修正的總體有限元方程組，是含所有待定未知量的封閉 方程組，採用適當的數值計算方法求解，可求得各節點的函數值。

有限體積法（Finite Volume Method）又稱為控制體積法。其基本思路是：將計算區域劃分為一系列不重復的控制體積，並使每個網格點周圍有一個控制體積；將待解的微分方程對每一個控制體積積分，便得出一組離散方程。其中的未知數是網格點上的因變數的數值。為了求出控制體積的積分，必須假定值在網格點之間的變化規律，即假設值的分段的分布的分布剖面。從積分區域的選取方法看來，有限體積法屬於加權剩餘法中的子區域法；從未知解的近似方法看來，有限體積法屬於採用局部近似的離散方法。簡言之，子區域法屬於有限體積發的基本方法。

有限體積法的基本思路易於理解，並能得出直接的物理解釋。離散方程的物理意義，就 是因變數在有限大小的控制體積中的守恆原理，如同微分方程表示因變數在無限小的控 制體積中的守恆原理一樣。 限體積法得出的離散方程，要求因變數的積分守恆對任意一組控制體積都得到滿足，對整個計算區域，自然也得到滿足。這是有限體積法吸引人的優點。有一些離散方法，例如有限差分法，僅當網格極其細密時，離散方程才滿足積分守恆；而有限體積法即使在粗網格情況下，也顯示出准確的積分守恆。就離散方法而言，有限體積法可視作有限單元法和有限差分法的中間物。有限單元法必須假定值在網格點之間的變化規律（既插值函數），並將其作為近似解。有限差分法只考慮網格點上的數值而不考慮值在網格點之間如何變化。有限體積法只尋求的結點值，這與有限差分法相類似；但有限體積法在尋求控制體積的積分時，必須假定值在網格點之間的分布，這又與有限單元法相類似。在有限體積法中，插值函數只用於計算控制 體積的積分，得出離散方程之後，便可忘掉插值函數；如果需要的話，可以對微分方程 中不同的項採取不同的插值函數。