有限元方法基本原理

一、数学基础

有限元方法最基本的数学工具是微积分和线性代数，因此需要在这些基础知识上做一些阐述。

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;
    }
}