简单有限元的C++编程

---

前言

接前文有限元编程示例matlab + C++，为了方便查看，将C++程序单独拿出来。具体问题分析可以查看前文

一、三连杆问题

问题描述:

C++程序:

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

MatrixXf Bar1D2Node_Stiffness(float E, float A, float L)
{
	//计算单元刚度矩阵，输入弹性模量E，横截面积A和长度L
	//输出单元刚度矩阵k(2X2) 

	MatrixXf k(2, 2);  //单元刚度矩阵的大小
	k << E * A / L, -E * A / L, -E * A / L, E* A / L;
	return k;
}

MatrixXf Bar1D2Node_Assembly(MatrixXf KK, MatrixXf k, int i, float j)
{
	//单元刚度矩阵的组装
	//输入单元刚度矩阵k，单元的节点编号i、j
	//输出整体刚度矩阵KK

	Vector2i Dof; //定义一个2*1的向量
	Dof << i, j;

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

int main()
{
	//弹性模量
	float E1 = 2e5;   float E2 = E1;      float E3 = E1;
	//面积
	float A3 = 0.06;  float A2 = A3 / 2;  float A1 = A3 / 3;
	//长度
	float L1 = 0.1;   float L2 = L1;      float L3 = L1;

	//分别计算三个单元的单刚
	MatrixXf k1(2, 2); MatrixXf k2(2, 2); MatrixXf k3(2, 2);//将三个单刚初始化
	k1 = Bar1D2Node_Stiffness(E1, A1, L1);
	k2 = Bar1D2Node_Stiffness(E2, A2, L2);
	k3 = Bar1D2Node_Stiffness(E3, A3, L3);

	//组装刚度矩阵
	MatrixXf KK(4, 4);    //总刚的大小
	KK.setZero(4, 4);
	KK = Bar1D2Node_Assembly(KK, k1, 1, 2);
	KK = Bar1D2Node_Assembly(KK, k2, 2, 3);
	KK = Bar1D2Node_Assembly(KK, k3, 3, 4);

	//求解节点位移
	MatrixXf kk(3, 3);
	kk = KK.block<3, 3>(0, 0);    //%节点位移u4 = 0;应用化1置0

	Vector3f P(-100, 0, 50);      //节点力
	Vector3f x(0, 0, 0);

	x = kk.lu().solve(P);         //LU分解求解线性方程组kk * x = P

	system("pause");
	return 0;
}

二、平面3节点的三角形单元

问题描述:

C++程序:

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

MatrixXf Tri2D3Node_Stiffness(float E, float mu, float t, float xi, float yi, float xj, float yj, float xm, float ym, int Id)
{
	//计算单元刚度矩阵，输入弹性模量E，泊松比NU，厚度t
	// 输入三个节点i、j、m的坐标xi,yi,xj,yj,xm,ym
	//输入平面问题性质指示参数ID(1为平面应力问题，2为平面应变)
	//输出单元刚度矩阵k(6X6) 
	float A, bi, bj, bm, ci, cj, cm;
	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, int i, int j, int m)
{
	//单元刚度矩阵的组装
	//输入单元刚度矩阵k，单元的节点编号i、j
	//输出整体刚度矩阵KK

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

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

int main()
{
	//初始物理量
	float E = 1e7;             //弹性模量
	float mu = 1.0 / 3;       //注意： 不能写成float mu = 1 / 3，这种情况输出结果为0
	float t = 0.1;             //厚度 
	int Id = 1;               //输入平面问题性质指示参数ID(1为平面应力问题，2为平面应变)，会给出不同的弹性矩阵D      

	MatrixXf k1(6, 6), k2(6, 6);     //单元刚度矩阵6*6
	k1 = Tri2D3Node_Stiffness(E, mu, t, 2, 0, 0, 1, 0, 0, Id);
	k2 = Tri2D3Node_Stiffness(E, mu, t, 0, 1, 2, 0, 2, 1, Id);
	cout << "k1:" << endl << k1 << endl;
	cout << "k2:" << endl << k2 << endl;

	//组装刚度矩阵
	MatrixXf KK(8, 8);    //总刚的大小
	KK.setZero(8, 8);     //0矩阵
	KK = Tri2D3Node_Assembly(KK, k1, 2, 3, 4);
	KK = Tri2D3Node_Assembly(KK, k2, 3, 2, 1);
	cout << "总体刚度矩阵 KK = " << endl << KK << endl;
	//计算节点向量   应用化1置0方法
	MatrixXf kk(4, 4);
	kk = KK.block<4, 4>(0, 0);    //节点位移u3,v3,u4,v4 = 0;应用化1置0可以把总刚矩阵简化为kk进行计算
	cout << "kk:" << endl << kk << endl;

	Vector4f p(0, -5000, 0, -5000);      //节点力
	Vector4f u(0, 0, 0, 0);
	u = kk.lu().solve(p);         //LU分解求解线性方程组kk * x = P
	cout << "节点位移 u = " << endl << u << endl;

	//支反力的计算
	VectorXf U(8);
	U << u, 0, 0, 0, 0;          //总体节点位移
	VectorXf P(8);
	P = KK * U;
	cout << "节点位移 U = " << endl << U << endl;
	cout << "支反力 P = " << endl << P << endl;

	system("pause");
	return 0;
}

总结

一维数组名称的用途：

二维数组定义的四种方式: