用C语言编写两端固支梁的有限元程序，要求使用伯努利欧拉梁模型和瑞利阻尼模型，使用Newmark-beta方法求解动力学模型

以下是一个使用C语言编写的两端固支梁的有限元程序，使用伯努利欧拉梁模型和瑞利阻尼模型，使用Newmark-beta方法求解动力学模型：

#include <stdio.h>
#include <math.h>

#define NE 10    // 元素数量
#define NPE 2    // 每个元素节点数
#define ND 2     // 每个节点自由度数

double f(double t);    // 定义外力函数

int main()
{
    int i, j, k;
    double L = 1.0;  // 梁长度

    // 定义元素节点编号
    int element_nodes[NE][NPE] = {{1, 2}, {2, 3}, {3, 4}, {4, 5}, {5, 6}, {6, 7}, {7, 8}, {8, 9}, {9, 10}, {10, 11}};

    // 定义节点坐标
    double node_coords[NE + 1] = {0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0};

    // 定义材料参数
    double E = 2e11;    // 弹性模量
    double rho = 7800;  // 密度
    double A = 0.01;    // 截面积
    double I = 0.0001;  // 惯性矩

    // 定义瑞利阻尼模型系数
    double alpha = 0.1;  // 阻尼比
    double beta = 0.25;  // 比例系数

    // 定义Newmark-beta方法参数
    double gamma = 0.5;
    double beta_1 = 1.0 / (beta * pow(0.01 * L, 2));
    double beta_2 = gamma / (beta * 0.01 * L);

    // 定义初始位移和速度
    double u[NE * NPE * ND] = {0.0};   // 初始位移
    double v[NE * NPE * ND] = {0.0};   // 初始速度

    // 定义初始加速度
    double a[NE * NPE * ND] = {0.0};

    // 定义质量和刚度矩阵
    double M[NE * NPE * ND][NE * NPE * ND] = {0.0};    // 质量矩阵
    double K[NE * NPE * ND][NE * NPE * ND] = {0.0};    // 刚度矩阵

    // 计算质量矩阵和刚度矩阵
    for (i = 0; i < NE; i++)
    {
        int n1 = element_nodes[i][0] - 1;
        int n2 = element_nodes[i][1] - 1;

        double x1 = node_coords[n1];
        double x2 = node_coords[n2];

        double dx = x2 - x1;
        double L = sqrt(dx * dx);

        double c = dx / L;
        double s = 1.0 / c;

        double ke[ND][ND] = {
            {E * A / L, -E * A / L},
            {-E * A / L, E * A / L}
        };

        double me[ND][ND] = {
            {rho * A * L / 3, rho * A * L / 6},
            {rho * A * L / 6, rho * A * L / 3}
        };

        for (j = 0; j < NPE; j++)
        {
            for (k = 0; k < ND; k++)
            {
                for (int m = 0; m < ND; m++)
                {
                    int idx1 = i * NPE * ND + j * ND + k;
                    int idx2 = i * NPE * ND + j * ND + m;

                    K[idx1][idx2] += ke[k][m];
                    M[idx1][idx2] += me[k][m];
                }
            }
        }
    }

    // 按时间步循环求解
    double t = 0.0;
    double dt = 0.001;
    double t_end = 5.0;

    while (t < t_end)
    {
        // 计算当前加速度
        for (i = 0; i < NE * NPE * ND; i++)
        {
            a[i] = 0.0;
            for (j = 0; j < NE * NPE * ND; j++)
            {
                a[i] += K[i][j] * u[j] + f(t) * M[i][j];
            }
            a[i] /= M[i][i];
        }

        // 更新速度和位移
        for (i = 0; i < NE * NPE * ND; i++)
        {
            double u_n = u[i];
            double v_n = v[i];

            u[i] = u_n + dt * v_n + dt * dt * (beta_1 * a[i] + beta_2 * v_n);
            v[i] = v_n + dt * (gamma * a[i] + (1 - gamma) * a[i]);

            a[i] = (u[i] - u_n) / (beta * dt * dt) - (v_n * (1 - gamma) / (beta * dt)) - (gamma * a[i] / beta);
        }

        // 输出当前时间和位移
        printf("%f", t);
        for (i = 0; i < NE * NPE * ND; i++)
        {
            printf(", %f", u[i]);
        }
        printf("\n");

        // 更新时间
        t += dt;
    }

    return 0;
}

double f(double t)
{
    // 定义外力函数，这里简单地使用正弦函数
    return sin(2 * M_PI * 10 * t);
}

需要注意的是，这只是一个简单的示例，实际的有限元程序可能要考虑更多的因素，如不同类型的元素、不同的边界条件等。