平面三节点有限元分析C++程序实现

---

前言

只是为方便学习，不做其他用途，matlab程序参考平面三角形单元有限元实现——有限元实践笔记（5）

一、有限元模型

相应的理论分析可参考平面三角形单元有限元实现——有限元实践笔记（5）
或者参考有限元分析的书三角形单元部分

二、C++程序实现

#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;

MatrixXf ele_xy(int i, MatrixXf ele, MatrixXf Node)
{
	//将第i+1个单元的节点对应的坐标  用3*2的矩阵存储起来
	//提取第i+1个单元的节点编号，并且将对应的三个节点编号对应的坐标用3*2的矩阵存储起来
	//ele:单元信息   Node：节点信息
	MatrixXf _ele(1, 3), b(1, 3), c(1, 3), d(1, 3);
	//matrix.block(i,j, p, q) : 表示返回从矩阵(i, j)开始，每行取p个元素，每列取q个元素所组成的临时新矩阵对象
	_ele = ele.block(i, 1, 1, 3);             //提取ele中第i+1个单元的3个节点编号   _ele：1行3列的向量
	b = Node.block(_ele(0, 0) - 1, 1, 1, 2);  //提取Node中第一个节点坐标
	c = Node.block(_ele(0, 1) - 1, 1, 1, 2);  //提取Node中第二个节点坐标
	d = Node.block(_ele(0, 2) - 1, 1, 1, 2);  //提取Node中第三个节点坐标
	MatrixXf node(3, 2);
	node << b, c, d;

	return node;
}

MatrixXf Tri2D3Node_Stiffness(float E, float mu, float t, MatrixXf ele_xy, int Id)
{
	//计算单元刚度矩阵，输入弹性模量E，泊松比NU，厚度t
	//矩阵ele_xy  对应单元的三个节点i、j、m的坐标xi,yi,xj,yj,xm,ym
	//输入平面问题性质指示参数ID(1为平面应力问题，2为平面应变)
	//输出单元刚度矩阵k(6X6) 
	float A, bi, bj, bm, ci, cj, cm, xi, xj, yi, yj, xm, ym;
	xi = ele_xy(0,0);
	yi = ele_xy(0, 1);
	xj = ele_xy(1, 0);
	yj = ele_xy(1, 1);
	xm = ele_xy(2, 0);
	ym = ele_xy(2, 1);
	A = (xi * (yj - ym) + xj * (ym - yi) + xm * (yi - yj)) / 2;
	bi = yj - ym;
	bj = ym - yi;
	bm = yi - yj;
	ci = xm - xj;
	cj = xi - xm;
	cm = xj - xi;
	MatrixXf B(3, 6);
	B << bi, 0, bj, 0, bm, 0,
		0, ci, 0, cj, 0, cm,
		ci, bi, cj, bj, cm, bm;
	B = B / (2 * A);

	MatrixXf D(3, 3);     // 弹性矩阵
	if (Id == 1)
	{
		D << 1, mu, 0,
			mu, 1, 0,
			0, 0, (1 - mu) / 2;
		D = (E / (1 - mu * mu)) * D;
	}
	else if (Id == 2)
	{
		D << 1 - mu, mu, 0,
			mu, 1 - mu, 0,
			0, 0, (1 - 2 * mu) / 2;
		D = (E / (1 + mu) / (1 - 2 * mu)) * D;
	}
	//B  矩阵转置：B.transpose()
	MatrixXf k(6, 6);  //单元刚度矩阵的大小
	k = t * A * B.transpose() * D * B;
	return k;
}

MatrixXf Tri2D3Node_Assembly(MatrixXf KK, MatrixXf k, MatrixXf _ele)
{
	//单元刚度矩阵的组装
	//输入单元信息，得到单元的节点编号i、j、m
	//输出整体刚度矩阵KK
	int i, j, m;
	i = _ele(0, 0);
	j = _ele(0, 1);
	m = _ele(0, 2);

	VectorXi Dof(6); //定义一个6*1的向量
	Dof << 2 * i - 2, 2 * i - 1, 2 * j - 2, 2 * j - 1, 2 * m - 2, 2 * m - 1;

	for (int n1 = 0; n1 < 6; n1++)
	{
		for (int n2 = 0; n2 < 6; n2++)
		{
			KK(Dof(n1), Dof(n2)) = KK(Dof(n1), Dof(n2)) + k(n1, n2);
			//注意：c++中数组和矩阵都是从0开始编号的
		}
	}
	return KK;
}

VectorXf Solve_1_Model(MatrixXf K, VectorXf U, VectorXf P, int len_U)
{
	//置“1”法求解方程组
	//K：总刚   U：初始化的节点位移向量  
	for (int j = 0; j < len_U; j++)
	{
		if (U(j) == 0)
		{
			K.block(j, 0, 1, K.cols()).setZero();    //将总刚矩阵K的行 置“0”  K.cols()：K矩阵的列长
			K.block(0, j, K.rows(), 1).setZero();    //将总刚矩阵K的列 置“0”  K.rows()：K矩阵的行长
			K(j, j) = 1;
			P(j) = 0;     //节点力置“0”
		}
	}

	// 节点位移
	U = K.lu().solve(P);         //LU分解求解线性方程组K * U = P
	return U;
}

int main()
{
	//初始物理量
	float E = 2.1e7;          //弹性模量
	float mu = 0.3;        
	float t = 0.025;          //厚度 
	int Id = 1;               //输入平面问题性质指示参数ID(1为平面应力问题，2为平面应变)，会给出不同的弹性矩阵D      
	
	//节点信息，第一列为节点编号，2-4列分别为x，y，z坐标
	MatrixXf Node(15,3);
	Node << 1, 0, 0, 
	2,   0.125,   0, 
	3,   0.25,    0, 
	4,   0.375,   0, 
	5,   0.5,     0,   
	6,   0,       0.125, 
	7,   0.125,   0.125, 
	8,   0.25,    0.125,
	9,   0.375,   0.125, 
	10,  0.5,     0.125, 
	11,  0,       0.25, 
	12,  0.125,   0.25, 
	13,  0.25,    0.25, 
	14,  0.375,   0.25, 
	15,  0.5,     0.25;
	
	//单元信息，第一列为单元编号，后面各列为单元上的节点号码
	MatrixXf ele(16, 4);
	ele << 1, 1, 2, 7,
		2, 1, 7, 6,
		3, 2, 3, 8,
		4, 2, 8, 7,
		5, 3, 4, 9,
		6, 3, 9, 8,
		7, 4, 5, 10,
		8, 4, 10, 9,
		9, 6, 7, 12,
		10, 6, 12, 11,
		11, 7, 8, 13,
		12, 7, 13, 12,
		13, 8, 9, 14,
		14, 8, 14, 13,
		15, 9, 10, 15,
		16, 9, 15, 14;

	int dof = Node.rows();       //节点数  15
	int n_ele = ele.rows();      //单元数  16
	MatrixXf k(6, 6);            //单刚矩阵
	MatrixXf K(dof * 2, dof * 2);        //总刚的大小
	K.setZero(dof * 2, dof * 2);         //矩阵置0

	for(int i = 0; i < n_ele; i++)
	{
		MatrixXf e_xy(3, 2);          
		e_xy = ele_xy(i, ele, Node);         //第i+1个单元的节点编号  对应的坐标用3*2矩阵存储
		k = Tri2D3Node_Stiffness(E, mu, t, e_xy, Id);   //第i+1个单刚矩阵
		MatrixXf _ele(1, 3);
		_ele = ele.block(i, 1, 1, 3);        //提取ele中第i+1个单元的3个节点编号
		K = Tri2D3Node_Assembly(K, k, _ele);
	}
	cout << "总刚矩阵 K = " << endl << K.block(15, 15, 15, 15) << endl;

	int len_U = 2 * Node.rows();        //节点位移矩阵长度
	//节点力
	VectorXf P(len_U);   P.setZero();
	//5、10、15  节点横向力
	P(8) = 4.6875;            P(18) = 9.375;            P(28) = 4.6875;           
	//初始节点位移矩阵，将非0的节点位移先  置“1”   方便之后用置“1”法求解方程
	VectorXf U(len_U);          U.setOnes();  //将矩阵置1
	U(0) = 0;U(1) = 0;U(10) = 0;U(11) = 0;U(20) = 0;U(21) = 0;
	
	U = Solve_1_Model(K, U, P,len_U);
	cout << "节点位移：U = " << endl << U << endl;
	// 支反力
	VectorXf F(len_U);
	F = K * U;
	cout << endl << "--------分割线-------- \n "  << endl << endl;
	cout << "支反力F = \n" << F << endl;

	system("pause");
	return 0;
}

该程序运行前需要添加eigen库
具体添加eigen库可参考C++中Eigen库的使用