C语言三梁平面框架结构的有限元分析

作者: 一路向后 | 来源:发表于2020-09-01 21:23 被阅读0次

C语言三梁平面框架结构的有限元分析
18、梁
根据两点经纬度求一条直线上的距离点/等分点（C#实现）（有问题，
C语言系列教程完整版（持续更新）
c语言四杆桁架结构的有限元分析
C语言(三)
第一章预备知识
三天学会C语言编程｜下篇
C，C++、C#三者有什么关系？
第一章预备知识

1.程序源码

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <malloc.h>
#include <math.h>

#define PI 3.141592654

//高斯消元法
int Gauss_Solve(double *x, double **A, double *b, int n)
{
        double mik;
        double s;
        int i, j, k;

        //消元
        for(k=0; k<n-1; k++)
        {
                if(!A[k][k])
                        return -1;

                for(i=k+1; i<n; i++)
                {
                        mik = A[i][k] / A[k][k];

                        for(j=k; j<n; j++)
                        {
                                A[i][j] = A[i][j] - mik * A[k][j];
                        }

                        b[i] = b[i] - mik * b[k];
                }
        }

        //printf("消元后的矩阵A: \n");

        //for(i=0; i<n; i++)
                //printf("%.4lf\n", b[i]);

        //回代
        x[n-1] = b[n-1] / A[n-1][n-1];

        for(k=n-2; k>=0; k--)
        {
                s = b[k];

                for(j=k+1; j<n; j++)
                {
                        s = s - A[k][j] * x[j];
                }

                x[k] = s / A[k][k];
        }

        return 0;
}

//3个方形矩阵乘法
int Matrix_Mul(double **A, double **B, double **C, double **D, int a)
{
    double **E = (double **)malloc(a*sizeof(double *));
    double *rE = (double *)malloc(a*a*sizeof(double));
    double sum;
    int i, j, k;

    for(i=0; i<a; i++)
    {
        E[i] = rE+i*a;
    }

    for(i=0; i<a; i++)
    {
        for(j=0; j<a; j++)
        {
            sum = 0;

            for(k=0; k<a; k++)
            {
                sum = sum + A[i][k] * B[k][j];
            }

            E[i][j] = sum;
        }
    }

    for(i=0; i<a; i++)
    {
        for(j=0; j<a; j++)
        {
            sum = 0;

            for(k=0; k<a; k++)
            {
                sum = sum + E[i][k] * C[k][j];
            }

            D[i][j] = sum;
        }
    }

    free(rE);
    free(E);

    return 0;
}

//计算刚度矩阵
int Beam2D2Node_Stiffness(double **k, double E, double I, double A, double L, double alpha)
{
    double **D = (double **)malloc(6*sizeof(double *));
    double **B = (double **)malloc(6*sizeof(double *));
    double **C = (double **)malloc(6*sizeof(double *));
    double *rA = (double *)malloc(36*sizeof(double));
    double *rB = (double *)malloc(36*sizeof(double));
    double *rC = (double *)malloc(36*sizeof(double));
    int i, j;

    for(i=0; i<6; i++)
    {
        D[i] = rA+i*6;
        B[i] = rB+i*6;
        C[i] = rC+i*6;
    }

    memset(rA, 0x00, sizeof(rA));
    memset(rB, 0x00, sizeof(rB));

    B[0][0] = cos(alpha);
    B[0][1] = sin(alpha);
    B[1][0] = -sin(alpha);
    B[1][1] = cos(alpha);
    B[2][2] = 1;
    B[3][3] = cos(alpha);
    B[3][4] = sin(alpha);
    B[4][3] = -sin(alpha);
    B[4][4] = cos(alpha);
    B[5][5] = 1;

    for(i=0; i<6; i++)
    {
        for(j=0; j<6; j++)
        {
            D[j][i] = B[i][j];
        }
    }

    C[0][0] =  E * A / L;
    C[0][1] =  0;
    C[0][2] =  0;
    C[0][3] = -E * A / L;
    C[0][4] =  0;
    C[0][5] =  0;

    C[1][0] =  0;
    C[1][1] =  12 * E * I / (L*L*L);
    C[1][2] =  6  * E * I / (L*L);
    C[1][3] =  0;
    C[1][4] = -12 * E * I / (L*L*L);
    C[1][5] =  6  * E * I / (L*L);

    C[2][0] =  0;
    C[2][1] =  6  * E * I / (L*L);
    C[2][2] =  4  * E * I /  L;
    C[2][3] =  0;
    C[2][4] = -6  * E * I / (L*L);
    C[2][5] =  2  * E * I /  L;

    C[3][0] = -E * A / L;
    C[3][1] =  0;
    C[3][2] =  0;
    C[3][3] =  E * A / L;
    C[3][4] =  0;
    C[3][5] =  0;

    C[4][0] =  0;
    C[4][1] = -12 * E * I / (L*L*L);
    C[4][2] = -6  * E * I / (L*L);
    C[4][3] =  0;
    C[4][4] =  12 * E * I / (L*L*L);
    C[4][5] = -6  * E * I / (L*L);

    C[5][0] =  0;
    C[5][1] =  6  * E * I / (L*L);
    C[5][2] =  2  * E * I /  L;
    C[5][3] =  0;
    C[5][4] = -6  * E * I / (L*L);
    C[5][5] =  4  * E * I /  L;

    Matrix_Mul(D, C, B, k, 6);

    free(rA);
    free(rB);
    free(rC);
    free(D);
    free(B);
    free(C);

    return 0;
}

//单元刚度矩阵的组装
int Beam2D2Node_Assembly(double **KK, double **k, int i, int j)
{
    KK[3*i+0][3*i+0] += k[0][0];
    KK[3*i+0][3*i+1] += k[0][1];
    KK[3*i+0][3*i+2] += k[0][2];
    KK[3*i+1][3*i+0] += k[1][0];
    KK[3*i+1][3*i+1] += k[1][1];
    KK[3*i+1][3*i+2] += k[1][2];
    KK[3*i+2][3*i+0] += k[2][0];
    KK[3*i+2][3*i+1] += k[2][1];
    KK[3*i+2][3*i+2] += k[2][2];

    KK[3*i+0][3*j+0] += k[0][3];
    KK[3*i+0][3*j+1] += k[0][4];
    KK[3*i+0][3*j+2] += k[0][5];
    KK[3*i+1][3*j+0] += k[1][3];
    KK[3*i+1][3*j+1] += k[1][4];
    KK[3*i+1][3*j+2] += k[1][5];
    KK[3*i+2][3*j+0] += k[2][3];
    KK[3*i+2][3*j+1] += k[2][4];
    KK[3*i+2][3*j+2] += k[2][5];

    KK[3*j+0][3*i+0] += k[3][0];
    KK[3*j+0][3*i+1] += k[3][1];
    KK[3*j+0][3*i+2] += k[3][2];
    KK[3*j+1][3*i+0] += k[4][0];
    KK[3*j+1][3*i+1] += k[4][1];
    KK[3*j+1][3*i+2] += k[4][2];
    KK[3*j+2][3*i+0] += k[5][0];
    KK[3*j+2][3*i+1] += k[5][1];
    KK[3*j+2][3*i+2] += k[5][2];

    KK[3*j+0][3*j+0] += k[3][3];
    KK[3*j+0][3*j+1] += k[3][4];
    KK[3*j+0][3*j+2] += k[3][5];
    KK[3*j+1][3*j+0] += k[4][3];
    KK[3*j+1][3*j+1] += k[4][4];
    KK[3*j+1][3*j+2] += k[4][5];
    KK[3*j+2][3*j+0] += k[5][3];
    KK[3*j+2][3*j+1] += k[5][4];
    KK[3*j+2][3*j+2] += k[5][5];

    return 0;
}

//边界条件处理
int Beam2D2Node_Condition(double **SS, double *b, double **KK, int *u, double *r, int n)
{
        int i, j;
        int p = 0, q = 0;

        for(i=0; i<n; i++)
        {
                if(u[i] == 0)
                {
                        b[p] = r[i];
                        q = 0;

                        for(j=0; j<n; j++)
                        {
                                if(u[j] == 0)
                                {
                                        SS[p][q] = KK[i][j];

                                        q++;
                                }
                        }

                        p++;
                }
        }

        return q;
}

//合并位移
int Beam2D2Node_Combine(double *r, double *b, int *u, int n)
{
        int i, j = 0;

        for(i=0; i<n; i++)
        {
                if(u[i] == 0)
                {
                        r[i] = b[j];
                        j++;
                }
        }

        return j;
}

//计算力矢量
int Beam2D2Node_Force(double *forces, double **KK, double *u, int n)
{
    int i = 0;
    int j = 0;

    for(i=0; i<n; i++)
    {
        forces[i] = 0;

        for(j=0; j<n; j++)
        {
            forces[i] += (KK[i][j] * u[j]);
        }
    }

    return 0;
}

//获取单元位移矢量
int Beam2D2Node_Offset(double *v, double *u, int i, int j, int k, int w, int s, int t)
{
        v[0] = u[i];
        v[1] = u[j];
        v[2] = u[k];
        v[3] = u[w];
        v[4] = u[s];
        v[5] = u[t];

        return 0;
}

//计算单元的几何矩阵
int Beam1D2Node_Strain(double *B, double x, double L, double y)
{
    double f = x / L;

    B[0] = -y * (12*f - 6) / (L*L);
    B[1] = -y * ( 6*f - 4) / L;
    B[2] =  y * (12*f + 6) / (L*L);
    B[3] = -y * ( 6*f - 2) / L;

    return 0;
}

//计算单元内某点的应力
int Beam2D2Node_Stress(double *stress, double E, double *B, double *u)
{
    *stress = E * (B[0] * u[0] + B[1] * u[1] + B[2] * u[2] + B[3] * u[3]);

    return 0;
}

//计算单元内某点的挠度
int Beam2D2Node_Deflection(double *v, double x, double L, double *u)
{
    double f = x / L;
    double N[4];

    N[0] = 1 - 3*f*f + 2*f*f*f;
    N[1] = L * (f - 2*f*f + f*f*f);
    N[2] = 3*f*f - 2*f*f*f;
    N[3] = L * (f*f*f - f*f);

    *v = u[0] * N[0] + u[1] * N[1] + u[2] * N[2] + u[3] * N[3];

    return 0;
}

//输出单元刚度矩阵
int Beam2Node_Printk(double **k)
{
        int i, j;

        printf("[");

        for(i=0; i<6; i++)
        {
                for(j=0; j<6; j++)
                {
                        if(j != 5)
                                printf("%.8lf ", k[i][j]);
                        else
                                printf("%.8lf", k[i][j]);
                }

                if(i != 5)
                        printf("\n");
        }

        printf("]\n");

        return 0;
}

//输出整体刚度矩阵
int Beam2Node_PrintK(double **k, int n)
{
        int i, j;

        printf("[");

        for(i=0; i<n; i++)
        {
                for(j=0; j<n; j++)
                {
                        if(j != n-1)
                                printf("%.8lf ", k[i][j]);
                        else
                                printf("%.8lf", k[i][j]);
                }

                if(i != n-1)
                        printf("\n");
        }

        printf("]\n");

        return 0;
}


int main()
{
    double rK[12][12];
    double rS[12][12];
    double *K[12];
    double *S[12];
    double b[12];
    double c[12];
    double r[12];
    double x[12];
    double rk[3][6][6];
    double *k[3][6];

    double E = 3e11;
    double I = 6.5e-7;
    double A = 6.8e-4;
    double L1 = 1.44;
    double L2 = 0.96;

    int u[12];
    int n;
    int i, j;

    memset(rK, 0x00, sizeof(rK));
    memset(rS, 0x00, sizeof(rS));
    memset(b, 0x00, sizeof(b));
    memset(c, 0x00, sizeof(c));
    memset(r, 0x00, sizeof(r));
    memset(x, 0x00, sizeof(x));
    memset(u, 0x00, sizeof(u));

        for(i=0; i<12; i++)
        {
                K[i] = rK[i];
                S[i] = rS[i];
        }

        for(i=0; i<3; i++)
        {
                for(j=0; j<6; j++)
                {
                        k[i][j] = rk[i][j];
                }
    }

    Beam2D2Node_Stiffness(k[0], E, I, A, L1, 0.00);
    Beam2D2Node_Stiffness(k[1], E, I, A, L2, PI/2);
    Beam2D2Node_Stiffness(k[2], E, I, A, L2, PI/2);

    Beam2D2Node_Assembly(K, k[0], 0, 1);
    Beam2D2Node_Assembly(K, k[1], 2, 0);
    Beam2D2Node_Assembly(K, k[2], 3, 1);

    //Beam2Node_Printk(k[0]);
    //Beam2Node_Printk(k[1]);
    //Beam2Node_Printk(k[2]);
    Beam2Node_PrintK(K, 12);

    u[6] = 1;
        u[7] = 1;
        u[8] = 1;
        u[9] = 1;
        u[10] = 1;
        u[11] = 1;

        r[0] = 3000;
        r[1] = -3000;
        r[2] = -720;
        r[3] = 0;
        r[4] = -3000;
        r[5] = 720;

        //加入边界条件
    n = Beam2D2Node_Condition(S, b, K, u, r, 12);

    Gauss_Solve(x, S, b, n);

    //Beam2Node_PrintK(S, n);

        //合并计算结果得出位移
    Beam2D2Node_Combine(c, x, u, 12);

    for(i=0; i<12; i++)
    {
        printf("%.8lf ", c[i]);
    }

    printf("\n");

    //计算各个分力
    Beam2D2Node_Force(r, K, c, 12);

    for(i=0; i<12; i++)
    {
        printf("%.8lf ", r[i]);
    }

    printf("\n");

    return 0;
}

2.编译源码

$ gcc -o Beam2D2Node Beam2D2Node.c

3.运行程序

./Beam2D2Node
[144311523.43750003 -0.04304201 1269531.25000000 -141666666.66666669 0.00000000 0.00000000 -2644856.77083333 0.04304201 1269531.25000000 0.00000000 0.00000000 0.00000000
-0.04304201 213283661.26543209 564236.11137150 0.00000000 -783661.26543210 564236.11111111 0.04304201 -212500000.00000000 0.00026039 0.00000000 0.00000000 0.00000000
1269531.25000000 564236.11137150 1354166.66666667 0.00000000 -564236.11111111 270833.33333333 -1269531.25000000 -0.00026039 406250.00000000 0.00000000 0.00000000 0.00000000
-141666666.66666669 0.00000000 0.00000000 144311523.43750003 -0.04304201 1269531.25000000 0.00000000 0.00000000 0.00000000 -2644856.77083333 0.04304201 1269531.25000000
0.00000000 -783661.26543210 -564236.11111111 -0.04304201 213283661.26543209 -564236.11085073 0.00000000 0.00000000 0.00000000 0.04304201 -212500000.00000000 0.00026039
0.00000000 564236.11111111 270833.33333333 1269531.25000000 -564236.11085073 1354166.66666667 0.00000000 0.00000000 0.00000000 -1269531.25000000 -0.00026039 406250.00000000
-2644856.77083333 0.04304201 -1269531.25000000 0.00000000 0.00000000 0.00000000 2644856.77083333 -0.04304201 -1269531.25000000 0.00000000 0.00000000 0.00000000
0.04304201 -212500000.00000000 -0.00026039 0.00000000 0.00000000 0.00000000 -0.04304201 212500000.00000000 -0.00026039 0.00000000 0.00000000 0.00000000
1269531.25000000 0.00026039 406250.00000000 0.00000000 0.00000000 0.00000000 -1269531.25000000 -0.00026039 812500.00000000 0.00000000 0.00000000 0.00000000
0.00000000 0.00000000 0.00000000 -2644856.77083333 0.04304201 -1269531.25000000 0.00000000 0.00000000 0.00000000 2644856.77083333 -0.04304201 -1269531.25000000
0.00000000 0.00000000 0.00000000 0.04304201 -212500000.00000000 -0.00026039 0.00000000 0.00000000 0.00000000 -0.04304201 212500000.00000000 -0.00026039
0.00000000 0.00000000 0.00000000 1269531.25000000 0.00026039 406250.00000000 0.00000000 0.00000000 0.00000000 -1269531.25000000 -0.00026039 812500.00000000]
0.00091766 -0.00001036 -0.00138737 0.00090119 -0.00001788 -0.00003883 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 
3000.00000000 -3000.00000000 -720.00000000 0.00000000 -3000.00000000 720.00000000 -665.78287259 2201.17836376 601.38524835 -2334.21712741 3798.82163624 1128.31159429

C语言三梁平面框架结构的有限元分析
1.程序源码 2.编译源码 3.运行程序
18、梁
1、梁的平面注写 2、梁的截面注写
根据两点经纬度求一条直线上的距离点/等分点（C#实现）（有问题，
此方法借助北京54平面坐标系计算。由C的源码转换而来。 C源码地址：C语言WGS84坐标转北京54坐标（高斯投影）...
C语言系列教程完整版（持续更新）
下面是C语言所有完整的教程 C语言教程（一）：程序开发理论基础 C语言教程（二）：开发环境搭建 C语言教程（三）：...
c语言四杆桁架结构的有限元分析
1.程序源码 2.编译源码 3.运行程序
C语言(三)
动态内存申请数组的长度是预先定义好的，在整个程序中固定不变，但是在实际编程中，所需的内存取决于实际输入的数据，无...
第一章预备知识
（简述） c++融合了三种不同的编程方式：c语言代表的过程性语言、c++在c语言基础上添加的类代表的面向对象的语言...
三天学会C语言编程｜下篇
本文是《三天学会C语言编程》的下篇，承接上一篇《三天学会C语言编程｜上篇》。本文重点介绍C语言的数据类型、表达式、...
C，C++、C#三者有什么关系？
C/C++/C#准确的说，是三种语言。 C语言既具有高级语言的特点，又具有汇编语言的特点。它可以作为工作系统设计语...
第一章预备知识
1.1 C++ 简介 C++融合了三种不同的编程方式：C 语言代表的过程性语言、添加到面向对象语言、C++ 模板支...