有限元编程作业——c语言三角形常应变单元计算程序

最新推荐文章于 2024-02-04 09:53:00 发布
最新推荐文章于 2024-02-04 09:53:00 发布
阅读量941 收藏 9
点赞数 5
文章标签： c语言 数学建模 算法
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/qwersdf_we/article/details/131332334
版权
该文章展示了一个C语言编写的程序，用于计算具有给定弹性模量、泊松比、节点坐标、荷载和边界条件的三角形单元的应变和应力。程序通过构建几何矩阵、应力弹性矩阵和单元刚度矩阵，实现了与有限元软件相当的计算结果。
摘要由CSDN通过智能技术生成

计算案例：

 

其中弹性模量设置为1，泊松比设置为0.2

节点坐标：1（3,0）2（1.5,0）3（0,0）4（3,-0.4）5（1.5,-0.4）6（0,-0.4）

等效节点荷载：1（0,-1），其他节点为0

边界条件：节点3,4的位移为0

C语言计算结果窗口：

 

 

 根据窗口提示输入初始条件，程序可实现任意节点数三角形单元计算，与有限元软件计算结果一致，下附C语言程序代码：

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

void print1(double **arr, int row, int col){  //输出矩阵

	int i = 0;
	int j = 0;
	for (i = 0; i < row; i++)
	{
		for (j = 0; j < col; j++)
		{
			printf("%6.2f", arr[i][j]);
		}
		printf("\n");
	} }

double **str (double x[3][2],int t,int E,double u,double*U,int m){
double x1, y1;
x1=x[0][0],y1=x[0][1];
double x2, y2;
x2=x[1][0],y2=x[1][1];
double x3, y3;
x3=x[2][0],y3=x[2][1];

 
double a1 = x2*y3 - x3*y2;
double a2 = x3*y1 - x1*y3;
double a3 = x1*y2 - x2*y1;
 
double b1 = y2 - y3;
double b2 = y3 - y1;
double b3 = y1 - y2;
 
double c1 = x3 - x2;
double c2 = x1 - x3;
double c3 = x2 - x1;
 
//算出该三角形的边长length1、2、3；
double length1 = sqrt((double)(x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2));
double length2 = sqrt((double)(x2 - x3)*(x2 - x3) + (y2 - y3)*(y2 - y3));
double length3 = sqrt((double)(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];
}
}
double strain[3];
memset(strain, 0.0, sizeof(strain));
for (int j = 0; j <3; j++){
for (int i = 0; i < 6; i++){
strain[j] +=B[j][i]*U[i];
}
}
double stress[3];
memset(stress, 0.0, sizeof(stress));
for (int j = 0; j <3; j++){
for (int i = 0; i < 6; i++){
stress[j] +=S[j][i]*U[i];
}
}printf("第%d个单元应变：\n",m);
	for (int j = 0; j < 3; j++)
		printf("%6.2f", strain[j]);
		printf("\n");
printf("第%d个单元应力：\n",m);
	for (int j = 0; j < 3; j++)
		printf("%6.2f", stress[j]);
		printf("\n");
return 0;}


double **gd (double x[3][2],int t,int E,double u){
double x1, y1;
x1=x[0][0],y1=x[0][1];
double x2, y2;
x2=x[1][0],y2=x[1][1];
double x3, y3;
x3=x[2][0],y3=x[2][1];

 
double a1 = x2*y3 - x3*y2;
double a2 = x3*y1 - x1*y3;
double a3 = x1*y2 - x2*y1;
 
double b1 = y2 - y3;
double b2 = y3 - y1;
double b3 = y1 - y2;
 
double c1 = x3 - x2;
double c2 = x1 - x3;
double c3 = x2 - x1;
 
//算出该三角形的边长length1、2、3；
double length1 = sqrt((double)(x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2));
double length2 = sqrt((double)(x2 - x3)*(x2 - x3) + (y2 - y3)*(y2 - y3));
double length3 = sqrt((double)(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;
 B1S = (double **)malloc(sizeof(double *) * 6);					
		for (int i = 0; i < 6; i++) {
 		B1S[i] = (double*)malloc(sizeof(double) *6);				
	}
for(int i=0;i<6;i++)
		for(int j=0;j<6;j++)
			B1S[i][j]=0;
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];
}}
for (int j = 0; j < 6; j++){
for (int i = 0; i < 6; i++){
B1S[i][j] *= (area*t);
}
}

return B1S;
}


int main() {
		//输入基本参数
int t=1;
int E;
printf("请输入弹性模量：\n");
scanf_s("%d", &E);

double u ;
printf("请输入泊松比：\n");
scanf_s("%lf", &u);

	int  i, j, line, column = 2;							//line表示节点数
	printf("请输入节点数：\n");
scanf_s("%d", &line);

	double **x;													//定义二位数组的首地址（首地址为二级指针）
 
	x = (double **)malloc(sizeof(double *) * line);					//给行指针分配内存空间	
	
	for (i = 0; i < line; i++) {
 
		x[i] = (double*)malloc(sizeof(double) * column);				//给每个行指针后面分配2个int的空间作为列
	}

	printf("请依次输入节点的x、y坐标:\n");
		for( i=0;i<line;i++)
		scanf_s("%lf%lf",&x[i][0],&x[i][1]);
			
//确定单元节点
int danyuan;
printf("请输入单元数：\n");
scanf_s("%d", &danyuan);
int **y;													
 
	y = (int **)malloc(sizeof(int *) * danyuan);					
	
	for (i = 0; i < line; i++) {
 
		y[i] = (int*)malloc(sizeof(int) * 3);				
	}

	printf("请依次输入组成单元的节点:\n");
		for( i=0;i<danyuan;i++)
		scanf_s("%d%d%d", &y[i][0],&y[i][1],&y[i][2]);
		
			
			//定义整体刚度矩阵	
	double **K;													
 
	K = (double **)malloc(sizeof(double *) * (2*line));						
	
	for (i = 0; i < (2*line); i++) {
 
		K[i] = (double*)malloc(sizeof(double) * (2*line));				
	}

	for(i=0;i<2*line;i++)
		for(j=0;j<2*line;j++)
			K[i][j]=0;

//计算刚度矩阵
	double p[3][2];
memset(p, 0.0, sizeof(p));

for(int m=0;m<danyuan;m++){
	for( i=0;i<3;i++){
	for( j=0;j<2;j++)
		p[i][j]=x[y[m][i]-1][j];}

double **ki;

	ki=gd ( p, t, E, u);

for( i=0;i<3;i++)
	for(j=0;j<3;j++){
K[2*y[m][i]-2][2*y[m][j]-2]+=ki[2*i][2*j];
K[2*y[m][i]-1][2*y[m][j]-2]+=ki[2*i+1][2*j];
K[2*y[m][i]-2][2*y[m][j]-1]+=ki[2*i][2*j+1];
K[2*y[m][i]-1][2*y[m][j]-1]+=ki[2*i+1][2*j+1];}
}
//输出刚度矩阵
printf("整体刚度矩阵：\n");
print1(K,2*line,2*line);

//输入等效荷载矩阵
double *P;													
	P = (double *)malloc(sizeof(double ) * (2*line));
	printf("请依次输入等效节点荷载:\n");
	for(i=0;i<line;i++){
				printf("节点%d:\n",i+1);
				scanf_s("%lf%lf", &P[2*i],&P[2*i+1]);}
//给定边界条件
double *U;													
	U = (double *)malloc(sizeof(double ) * (2*line));
	int RANK=2*line;
	printf("请输入边界条件(位移为0时输入0，不为0时输入1):\n");
	for(i=0;i<2*line;i++)
		scanf_s("%lf",&U[i]);
for(i=0;i<2*line;i++){
	if(U[i]==0){
for(int m=0;m<2*line;m++)
	{K[m][i]=0.0;
	K[i][m]=0.0;}
	P[i]=0.0;
K[i][i]=1.0;}
}

for(i=0;i<2*line;i++)	
	U[i]=0;

//计算节点位移
	double get_A = 0.0;
	printf("利用以上A与b组成的增广阵进行高斯消去法计算方程组\n");
	for (int i = 1; i < RANK; i++)
	{
		for (int j = i; j<RANK; j++)
		{
			get_A = K[j][i - 1] / K[i - 1][i - 1];
			P[j] = P[j] - get_A * P[i - 1];
			for (int k = i-1; k < RANK; k++)
			{
				K[j][k] = K[j][k] - get_A * K[i-1][k];
			}
		}
	}
	printf("顺序消元后的上三角系数增广矩阵如下\n");
	for (int i = 0; i < RANK; i++)
	{
		for (int j = 0; j<RANK; j++)
		{
			printf("%g\t", K[i][j]);
		}
		printf("    %g", P[i]);
		printf("\n");
	}
	printf("节点位移：\n");

	for (int i = 0; i < RANK; i++)
	{
		double get_x = 0.0;
		for (int j = 0; j < RANK; j++)
		{
			get_x = get_x + K[RANK-1-i][j]*U[j];
		}
		U[RANK - 1 - i] = (P[RANK - 1 - i] - get_x + K[RANK - 1 - i][RANK - 1 - i] * U[RANK - 1 - i]) / K[RANK - 1 - i][RANK - 1 - i];
	}
	for (int i = 0; i < RANK; i++)
	{
		printf("x%d = %g\n", i + 1, U[i]);
	}

//计算单元应变应力
for(int m=0;m<danyuan;m++){
	for( i=0;i<3;i++){
	for( j=0;j<2;j++)
		p[i][j]=x[y[m][i]-1][j];}
double *v;
v = (double *)malloc(sizeof(double ) * (6));
for(i=0;i<3;i++)
	{v[2*i]=U[2*y[m][i]-2];
	v[2*i+1]=U[2*y[m][i]-1];
}
str( p, t, E, u,v,m+1);
}
	system("pause");
return 0;
}

qwersdf_we
关注
  • 5
    点赞
  • 踩
  • 9
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
有限元之编程示例三：有限元求抛物型偏微分方程
欢迎各位大佬光临“寒舍”！请多多指点，共同进步！
06-05 290
主要对有限元法求解抛物型偏微分方程边值问题进行C++编程。
空间梁单元C语言有限元编程,C语言三梁平面框架结构的有限元分析
weixin_30430517的博客
05-18 531
1.程序源码#include #include #include #include #include #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{if(!A[k][k])retu...
c语言求解有限元中常应变三角形单元分析问题
鱼木戈大
07-24 5374
1,程序设计目的   巩固和加深对有限元基本理论的理解，掌握程序设计的基本方法,技能及程序语言的编写及编译，提升求解实际工程问题。 2,程序说明 运用有限元方法中三角形常应变单元解平面问题的计算主程序。其步骤如下： （1）由结点位移及坐标确定出单元内部各点位移。 （2）由单元内部各点的位移确定出对应的应变。 （3）根据应力——应变的关系由单元应变确定出相应的应力。 （4）根据虚功方程...
用C++写一个简单的有限元程序
weixin_42584586的博客
02-13 560
下面是一个用C语言写的简单的有限元程序的示例： #include <stdio.h> #include <math.h> // 定义一个有限元函数 double finite_element(double x) { return x * x; } int main() { // 定义一个变量，存储积分结果 double result = 0.0;...
有限元六面体单元 C语言源代码.rar
04-15
有限元大作业，用C语言编写有限六面体单元源代码，压缩包内包括有限元六面体单元C语言代码的主程序以及计算行列式代码！
平面三节点有限元分析C++程序实现
mw_1422102031的博客
11-11 1255
平面三角形单元有限元实现 C++
有限元一维算例c语言,有限元方法入门（1）：一维算例
weixin_29251131的博客
05-24 1543
8种机械键盘轴体对比本人程序员，要买一个写代码的键盘，请问红轴和茶轴怎么选？有限元方法多用于二维、三维的区域上，相比于有限差分方法，对非规则网格的适应性较好。而在一维区域上，有限元方法与有限差分方法相比，没有太大的优势。但从一维的例子讲起，有助于由浅及深的了解有限元的思想。有限元方法多用于椭圆方程的求解，考虑最为常见的椭圆型方程——Poisson方程，如下：其中 $Omega$ 是 $x$ 区域...
爆炸三角形如何用计算机实现,爆炸三角形原理
weixin_33308985的博客
07-01 605
原标题：爆炸三角形原理国家石油化学工业的高速发展为企业带来巨大经济效益与社会效益的同时，各种安全隐患及事故也随之不断产生，给社会带来了极大的危害。作为保护生产安全的设备——— 安全栅 ，以前只在重要的场合应用, 现在已经作为设计的标准规范进行推广。1 安全栅的理论基础众所周知,化工企业属于高危企业，在生产运行中易发生爆炸等重大事故。在现场同时满足下列 3 个条件时，爆炸就会发生：①易爆物质如易爆气...
CST常应变三角形单元及其Matlab代码
qq_39642029的博客
05-08 1131
function k=Triangle2DD3Node_Stiffness(E,NU,t,xi,yi,xj,yj,xm,ym,ID) %该函数计算单元的刚度矩阵 %输入弹性模量E、泊松比NU和厚度t %输入3个节点i,j,m的坐标xi,yi,xj,yj,xm,ym %输入平面问题性质指示参数ID(1为平面应力，2为平面应变) %输出单元刚度矩阵k(6x6) %-------------------------------------------------- A=(xi*(yj-ym)+x...
有限元-二维有限元编程（矩形区域、三角剖分）
栀彧
07-13 5108
有限元、矩形区域、三角剖分
基于MATLAB编程实现的平面三角形单元通用有限元分析程序
Hardly_coding的博客
01-11 1925
基于MATLAB语言实现的有限元求解程序，从INP文件中读取数据信息，使用平面三角形单元划分网格，求解结果和ABAQUS大型商业软件的分析结果相比精度很好
有限元四面体单元 C语言源代码.rar
04-15
有限元大作业，用C语言编写有限四面体单元源代码，压缩包内包括有限元四面体单元C语言代码的主程序以及计算行列式代码！
C语言与有限元程序设计
12-25
本文通过与民件以入语言的比较, 说明利用语言开发有限元程序的可行性。 结果表明, 语言不但完全可以实现皿的语法功能, 拥有更强大的实用编程环境, 同时由于提供指针变量、动态内存分配函数和结构变量, 使其编写的软件在维护性、可读性和 内存利用效率方面具有明显的优势。文章最后还将语言与最近出现的’ 召作了简 单比较。
有限元的C语言实现计算
09-19
做的比较简单的C实现有限元的计算，由于非常简单就不收资源，大家学习探讨用，
３维矢量有限元C程序
04-06
３维矢量有限元C语言程序，针对矢量有限元方面的计算，为大家提供一个例子
三角形单元有限元程序
05-10
有限元三角形单元程序，很经典！希望大家能够互相学习与分享。
有限元编程学习方法及相关资料链接
最新发布
u010542847的博客
02-04 400
学习有限元编程是一个深入理解有限元方法、数值分析和结构力学的过程。以下是一个分阶段的学习方法，每个阶段提供相关资源链接和建议，以帮助你逐步掌握有限元编程的基础和高级知识。
有限元c语言编程大作业,C语言及有限元程序设计.pdf
weixin_42168830的博客
05-16 1224
C语言与有限元程序设计叶 又 戚 燕(上海交通大攀国家棋其C几D工程中心 上海加0030) (上海大学计算中心 上海匆印力摘 要 本文通过与FORTRAN77语官的比较，说明利用C语言开发有限元程序的可行性。结果表明,C语官不但完全可以实现FORTRAN77的语法功能，拥有更强大的实用编程环境，同时由于提供指针变f、动态内存分配函数和结构...
MATLAB实现的3节点三角形单元有限元程序
"3结点三角形单元有限元程序MATLAB语言是一个用于解决结构力学问题的计算程序，它利用有限元方法分析二维平面上的结构。程序包含6个主要部分，分别是主程序、总刚度计算、结点位移求解、线性方程组求解、应力应变...
评论 1
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值