| c语言求解有限元中常应变三角形单元分析问题 | 您所在的位置：网站首页 › 使用c语言求解实际问题的步骤是什么 › c语言求解有限元中常应变三角形单元分析问题 |
c语言求解有限元中常应变三角形单元分析问题
|
1,程序设计目的 巩固和加深对有限元基本理论的理解，掌握程序设计的基本方法,技能及程序语言的编写及编译，提升求解实际工程问题。 2,程序说明 运用有限元方法中三角形常应变单元解平面问题的计算主程序。其步骤如下： （1）由结点位移及坐标确定出单元内部各点位移。 （2）由单元内部各点的位移确定出对应的应变。 （3）根据应力——应变的关系由单元应变确定出相应的应力。 （4）根据虚功方程由应力确定出结点力。 （5）由单元的结点位移求出单元的结点力，据此求出单元的刚度矩阵。 3,应用实例 源程序： //求图3-6所示等腰直角三角形单元的刚度矩阵k,设厚度t=1,弹性模量E=1； #include #include #include #include #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; }运行上述代码； 输入三个节点的坐标，输入泊松比；
|
【本文地址】
公司简介
联系我们