c语言求解有限元中常应变三角形单元分析问题

1,程序设计目的

  巩固和加深对有限元基本理论的理解，掌握程序设计的基本方法,技能及程序语言的编写及编译，提升求解实际工程问题。

2,程序说明

运用有限元方法中三角形常应变单元解平面问题的计算主程序。其步骤如下：

（1）由结点位移及坐标确定出单元内部各点位移。

（2）由单元内部各点的位移确定出对应的应变。

（3）根据应力——应变的关系由单元应变确定出相应的应力。

（4）根据虚功方程由应力确定出结点力。

（5）由单元的结点位移求出单元的结点力，据此求出单元的刚度矩阵。

3,应用实例

源程序：

//求图3-6所示等腰直角三角形单元的刚度矩阵k,设厚度t=1,弹性模量E=1；

#include <stdio.h>
#include <Windows.h>
#include <math.h>
#include <stdlib.h>
#pragma warning(disable:4996)

menu(){
printf("-------------------------------------------------------------------------------\n");
printf("--------------------------Welcome to 有限元程序设计----------------------------\n");
printf("-------------------------------------------------------------------------------\n");
printf("-------------------------------------------------------------------------------\n");
printf("-------------------------------------------------------------------------------\n");
printf("-------------------该程序的目的为求解三角形单元的刚度矩阵Ke!-------------------\n");
printf("-------------------------------------------------------------------------------\n");
printf("-------------------------------------------------------------------------------\n");
printf("--------------------------假设厚度t=1，弹性模量E=1-----------------------------\n");
printf("-------------------------------------------------------------------------------\n");
printf("----------------------------按任意键进入程序设计！-----------------------------\n");
printf("-------------------------------------------------------------------------------\n");
printf("-------------------------------------------------------------------------------\n");
char start = 0;
scanf("%c", &start);
}
int main()
{
static int E = 1;
static int t = 1;
//先创建三个节点的x、y坐标

menu();
system("cls");
//结点1
int x1, y1;
//结点2
int x2, y2;
//结点3
int x3, y3;
printf("请输入第一个节点的x、y坐标:\n");
scanf("%d%d", &x1,&y1);

printf("请输入第二个节点的x、y坐标:\n");
scanf("%d%d", &x2,&y2);

printf("请输入第三个节点的x、y坐标:\n");
scanf("%d%d", &x3,&y3);
//输入泊松比；
double u = 0.0;
printf("请输入泊松比：\n");
scanf("%f", &u);

int a1 = x2*y3 - x3*y2;
int a2 = x3*y1 - x1*y3;
int a3 = x1*y2 - x2*y1;

int b1 = y2 - y3;
int b2 = y3 - y1;
int b3 = y1 - y2;

int c1 = x3 - x2;
int c2 = x1 - x3;
int c3 = x2 - x1;

//算出该三角形的边长length1、2、3；
double length1 = sqrt((x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2));
double length2 = sqrt((x2 - x3)*(x2 - x3) + (y2 - y3)*(y2 - y3));
double length3 = sqrt((x1 - x3)*(x1 - x3) + (y1 - y3)*(y1 - y3));
//算出该三角形的周长per
double per = (length1 + length2 + length3) / 2;
//算出该三角形的面积area；
double area = sqrt(per*(per - length1)*(per - length2)*(per - length3));

//建立几何矩阵（用二维数组实现）
double B[3][6];
memset(B, 0.0, sizeof(B));
//二维数组初始化
double cd = 1 / (2 * area);
//根据课本中的3-9公式
//第一行
B[0][0] = b1*cd;
B[0][2] = b2*cd;
B[0][4] = b3*cd;
//第二行
B[1][1] = c1*cd;
B[1][3] = c2*cd;
B[1][5] = c3*cd;
//第三行
B[2][0] = c1*cd;
B[2][1] = b1*cd;
B[2][2] = c2*cd;
B[2][3] = b2*cd;
B[2][4] = c3*cd;
B[2][5] = b3*cd;
//完成几何矩阵的创建；

//由于u=0，平面应力与平面问题的D相同；
//计算应力弹性矩阵
double D[3][3] = { { 1.0, u, 0.0 }, \
{u, 1.0, 0.0}, \
{0.0, 0.0, (1 - u) / 2}
};
for (int i = 0; i < (sizeof(D) / sizeof(D[0])); i++){
D[0][i] *= (E / (1 - u*u));
D[1][i] *= (E / (1 - u*u));
D[2][i] *= (E / (1 - u*u));
}//求得最终的应力弹性矩阵；

//计算应力转换矩阵
double S[3][6];
//数组初始化；
memset(S, 0.0, sizeof(S));

//S=D*B;
for (int j = 0; j < (sizeof(D) / sizeof(D[0])); j++){
for (int i = 0; i < (sizeof(D) / sizeof(D[0])); i++){
S[j][0] += D[j][i] * B[i][0];
S[j][1] += D[j][i] * B[i][1];
S[j][2] += D[j][i] * B[i][2];
S[j][3] += D[j][i] * B[i][3];
S[j][4] += D[j][i] * B[i][4];
S[j][5] += D[j][i] * B[i][5];
}
}

//将B、S及t，面积area带入公式得单元刚度矩阵Ke；

//求得B的转置矩阵B1；
double B1[6][3];
memset(B1, 0.0, sizeof(B1));
for (int i = 0; i < (sizeof(B1) / sizeof(B1[0])); i++){
B1[i][0] = B[0][i];
B1[i][1] = B[1][i];
B1[i][2] = B[2][i];
}

//求单元刚度矩阵；
double B1S[6][6];
memset(B1S, 0.0, sizeof(B1S));
for (int j = 0; j < 6; j++){
for (int i = 0;  i< 3; i++){
B1S[j][0] += B1[j][i] * S[i][0];
B1S[j][1] += B1[j][i] * S[i][1];
B1S[j][2] += B1[j][i] * S[i][2];
B1S[j][3] += B1[j][i] * S[i][3];
B1S[j][4] += B1[j][i] * S[i][4];
B1S[j][5] += B1[j][i] * S[i][5];
}
}
printf("求得该三角形单位刚度矩阵Ke为：\n");
printf("-------------------------------------------------------------------------------\n");
for (int j = 0; j < (sizeof(B1S) / sizeof(B1S[0])); j++){
for (int i = 0; i < (sizeof(B1S) / sizeof(B1S[0])); i++){
printf("%6.2f ",B1S[i][j] *= (area*t));//为最终结果！
}
printf("\n");
}
printf("-------------------------------------------------------------------------------\n");
system("pause");
return 0;
}

运行上述代码；

输入三个节点的坐标，输入泊松比；