有限元方法基本原理

一、數學基礎

有限元方法最基本的數學工具是微積分和線性代數，因此需要在這些基礎知識上做一些闡述。

1.微積分

微積分是求解物理問題的基礎，其中積分和微分是必不可少的工具。在有限元方法中，微積分主要用於對微分方程建立模型並推導其解析解或數值解。

<script type="math/tex">f(x) = \int_a^b F(x, t) dt</script>

在上面的代碼中，f(x)表示某個函數在x點處的值，F(x, t)是一個定義在區域[a,b]上的函數。有限元方法中，幾何區域一般使用三角形劃分來進行近似，因此需要將一維的積分轉化為二維的。

2.線性代數

在有限元方法中，線性代數主要用於解線性方程組，例如在有限元分析中需要求解如下的方程組：

<script type="math/tex">
KU = F
</script>

其中，K是一個$n \times n$的剛度矩陣，U是未知向量，F是力向量。解這個方程組可以得到結構的位移解。

二、數學模型

有限元方法中的數學模型一般來自物理系統的微分方程，將微分方程進行離散化，形成有限元方法求解的模型。

1.一維彈性力學模型

一維彈性力學模型是最基礎的有限元方法模型，對於一個線彈性材料，它的微分方程可以表示為：

<script type="math/tex">
\frac{d^2 u}{dx^2} + \frac{F}{AE} = 0
</script>

其中，u是位移，x是位置，F是力，A是截面面積，E是彈性模量。將其進行離散化，即可以得到相應的有限元方法模型。

2.二維熱傳導模型

二維熱傳導模型是用有限元方法求解熱傳導問題時採用的模型。它的微分方程可以表示為：

<script type="math/tex">
\frac{\partial u}{\partial t} - \alpha (\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}) = 0
</script>

其中，u是溫度，t是時間，$\alpha$是熱傳導係數。將其進行離散化，可以得到相應的有限元方法模型。

三、有限元離散化及求解

將微分方程離散化後，需要進行有限元離散化，建立節點、單元等概念，並將參考單元映射到實際單元上。常用的有限元類型有線性三角形單元、線性四邊形單元等。

1.線性三角形單元

線性三角形單元是最基礎的有限元方法單元，它由三個節點構成，具有良好的計算效率和精度，可用於求解大多數的應用問題。

<script type="math/tex">
\vec{u}(x) = N_1(x)\vec{u_1}+N_2(x)\vec{u_2}+N_3(x)\vec{u_3}
</script>

其中，$N_1(x),N_2(x),N_3(x)$是形函數，$\vec{u_1},\vec{u_2},\vec{u_3}$是節點上的位移向量。

2.求解過程

有了離散化後的模型和節點、單元數據，就可以開始進行有限元方法的求解了。整個過程可以分為如下步驟：

（1）剛度矩陣和力向量的組裝

需要對每個單元求其剛度矩陣和力向量，然後將其組裝成整個系統的剛度矩陣和力向量。

<script type="math/tex">
K_{ij} = \sum_{k=1}^{ne} k^{(k)}_{ij}, F_{i} = \sum_{k=1}^{ne} f^{(k)}_{i}
</script>

其中，K是整個系統的剛度矩陣，F是整個系統的力向量，$k^{(k)}$是單元k的剛度矩陣，$f^{(k)}$是單元k的力向量。

（2）加載和約束條件的處理

根據實際問題的邊界條件和約束條件，需要對載荷向量和剛度矩陣進行相應的調整。

<script type="math/tex">
K_{ij} = K_{ij} + K_{bc}, F_{i} = F_{i} - F_{bc}
</script>

其中，$K_{bc}$是邊界條件影響的剛度矩陣，$F_{bc}$是邊界條件影響的力向量。

（3）求解方程組

將上述處理後的剛度矩陣和力向量帶入方程組，使用高斯消元法、迭代法等求解數值解。

<script type="math/tex">
KU = F
</script>

四、代碼實現

// 計算單元剛度矩陣
void ElementStiffness(double k[3][3], double area, double lambda, double mu, double x1, double y1, double x2, double y2, double x3, double y3)
{
    double b[3], c[3], a[3];
    a[0] = y2 - y3;
    a[1] = y3 - y1;
    a[2] = y1 - y2;
    b[0] = x3 - x2;
    b[1] = x1 - x3;
    b[2] = x2 - x1;
    c[0] = x2 * y3 - x3 * y2;
    c[1] = x3 * y1 - x1 * y3;
    c[2] = x1 * y2 - x2 * y1;
    double db[2][2], dc[2][2];
    for (int i = 0; i < 2; i++)
        for (int j = 0; j < 2; j++)
        {
            db[i][j] = b[i + 1] * b[j + 1];
            dc[i][j] = c[i + 1] * c[j + 1];
        }
    for (int i = 0; i < 3; i++)
        for (int j = 0; j < 3; j++)
            k[i][j] = area * (lambda * (b[i] * b[j] + a[i] * a[j]) / (4 * area) + mu * (db[i % 2][j % 2] + dc[i % 2][j % 2]) / (2 * area));
}

// 組裝剛度矩陣和載荷向量
void AssembleStiffnessAndLoad(int n, double *xCoord, double *yCoord, double mu, double lambda, double *f, double *stiffMat, double *loadVec)
{
    double k[3][3];
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
        {
            stiffMat[i * n + j] = 0;
            loadVec[i] = 0;
        }
    for (int iE = 0; iE < n - 2; iE++)
    {
        double x1 = xCoord[iE], y1 = yCoord[iE];
        double x2 = xCoord[iE + 1], y2 = yCoord[iE + 1];
        double x3 = xCoord[iE + 2], y3 = yCoord[iE + 2];
        double area = fabs((x1 - x3) * (y2 - y3) - (x2 - x3) * (y1 - y3)) / 2;
        ElementStiffness(k, area, lambda, mu, x1, y1, x2, y2, x3, y3);
        stiffMat[iE * n + iE] += k[0][0];
        stiffMat[iE * n + iE + 1] += k[0][1];
        stiffMat[iE * n + iE + 2] += k[0][2];
        stiffMat[(iE + 1) * n + iE] += k[1][0];
        stiffMat[(iE + 1) * n + iE + 1] += k[1][1];
        stiffMat[(iE + 1) * n + iE + 2] += k[1][2];
        stiffMat[(iE + 2) * n + iE] += k[2][0];
        stiffMat[(iE + 2) * n + iE + 1] += k[2][1];
        stiffMat[(iE + 2) * n + iE + 2] += k[2][2];
        loadVec[iE] += f[iE] * area / 3;
        loadVec[iE + 1] += f[iE] * area / 3;
        loadVec[iE + 2] += f[iE] * area / 3;
    }
}