C++实现有限元二维杆单元计算 Bar2D2Node类(纯自研非套壳)

本系列文章致力于实现“手搓有限元，干翻Ansys的目标”，基本框架为前端显示使用QT实现交互，后端计算采用Visual Studio C++。

QT软件界面

具体软件操作可查看下方视频哦。也可以点击这里直接跳转。

直接干翻Ansys？小伙自研有限元

1、Bar2D2Node类

Bar2D2Node类用来实现有限元中二维杆单元的计算。整体架构如下：

Bar2D2Node类架构图

1.1、public function

1.1.1、构造函数与析构函数

构造函数用来初始化二维杆单元基本信息，包括ID、杆单元起始点、结束点、杨氏模量、横截面积的关键数值；析构函数用来释放内存。

Bar2D2Node.h函数声明文件

//***********************构造函数析构函数***********************//
/*
函数名称:		无参构造函数
*/
Bar2D2Node();/*
函数名称:		有参构造函数
id:				ID
*p0:			起始点
*p1:			结束点
E:				弹性模量
A:				横截面积
*/
Bar2D2Node(int id, Point2D* p0, Point2D* p1, double E, double A);/*
函数名称:		析构函数
*/
~Bar2D2Node();

Bar2D2Node.cpp函数实现文件

//无参构造函数
Bar2D2Node::Bar2D2Node()
{
}//有参构造函数
Bar2D2Node::Bar2D2Node(int id, Point2D* p1, Point2D* p2, double E, double A)
{this->ID = id;this->m_Point0 = p1;this->m_Point1 = p2;this->m_E = E;this->m_A = A;
}//析构函数
Bar2D2Node::~Bar2D2Node()
{
}

1.1.2、设置数值接口函数

设置数值接口函数可以单独设置二维杆单元的ID、杆单元起始点、结束点、杨氏模量、横截面积。

Bar2D2Node.h函数声明文件

//***********************设置数值接口函数***********************//
/*
函数名称:		设置Bar2D2Node杆单元
id:				ID
*p0:			起始点
*p1:			结束点
E:				弹性模量
A:				横截面积
*/
virtual void		SetBar(int id, Point2D* p0, Point2D* p1, double E, double /*
函数名称:		设置Bar2D2Node起始点与结束点
*p0:			起始点
*p1:			结束点
*/
virtual void		SetPoint(Point2D* p0, Point2D* p1);/*
函数名称:		设置Bar2D2Node ID
id:				ID
*/
void		SetID(int id);/*
函数名称:		设置Bar2D2Node弹性模量
E:				弹性模量
*/
void		SetE(double E);/*
函数名称:		设置Bar2D2Node截面面积
A:				横截面积
*/
void		SetA(double A);

Bar2D2Node.cpp函数实现文件

//设置Bar2D2Node杆单元
void Bar2D2Node::SetBar(int id, Point2D* p0, Point2D* p1, double E, double A)
{this->ID = id;this->m_Point0 = p0;this->m_Point1 = p1;this->m_E = E;this->m_A = A;
}//设置Bar2D2Node起始点与结束点
void Bar2D2Node::SetPoint(Point2D* p1, Point2D* p2)
{this->m_Point0 = p1;this->m_Point1 = p2;
}//设置Bar2D2Node ID
void Bar2D2Node::SetID(int id)
{this->ID = id;
}//设置Bar2D2Node弹性模量
void Bar2D2Node::SetE(double E)
{this->m_E = E;
}//设置Bar2D2Node截面面积
void Bar2D2Node::SetA(double A)
{this->m_A = A;
}

1.1.3、获取数值接口函数

获取数值接口函数不仅可以获取到二维杆单元的ID、杆单元起始点、结束点、杨氏模量、横截面积这些初始化变量。还可以获取应力、应变这些有限元计算结果变量，有的同学可能会比较疑惑，最基本的位移变量怎么没有获取呢？这是因为在二维杆单元的起始点、结束点（Point2D类）中已经包含节点位置信息。点击此处可以查看Point2D类的相关介绍。

Bar2D2Node.h函数声明文件

//***********************获取数值接口函数***********************//
/*
函数名称:		获取Bar2D2Node起始点
*/
Point2D*	GetPoint0();/*
函数名称:		获取Bar2D2Node终止点
*/
Point2D*	GetPoint1();/*
函数名称:		获取Bar2D2Node单元ID
*/
int			GetID();/*
函数名称:		获取Bar2D2Node单元弹性模量
*/
double		GetE();/*
函数名称:		获取Bar2D2Node单元面积
*/
double		GetA();/*
函数名称:		获取Bar2D2Node单元刚度矩阵
*/
Matrix		GetK();/*
函数名称:		获取Bar2D2Node单元应变
*/
double		GetEpsilon();/*
函数名称:		获取Bar2D2Node单元应力
*/
double		GetSigama();

Bar2D2Node.cpp函数实现文件

//获取Bar2D2Node起始点
Point2D* Bar2D2Node::GetPoint0()
{return this->m_Point0;
}//获取Bar2D2Node终止点
Point2D* Bar2D2Node::GetPoint1()
{return this->m_Point1;
}//获取Bar2D2Node单元ID
int Bar2D2Node::GetID()
{return this->ID;
}//获取Bar2D2Node单元弹性模量
double Bar2D2Node::GetE()
{return this->m_E;
}//获取Bar2D2Node单元面积
double Bar2D2Node::GetA()
{return this->m_A;
}//获取Bar2D2Node单元刚度矩阵
Matrix Bar2D2Node::GetK()
{return this->m_MatK;
}//获取Bar2D2Node单元应变
double Bar2D2Node::GetEpsilon()
{return this->m_Epsilon;
}//获取Bar2D2Node单元应力
double Bar2D2Node::GetSigama()
{return this->m_Sigama;
}

1.1.4、计算函数

此部分函数为有限元计算核心，可以计算二维杆单元全局刚度矩阵、应变、应力关键力学信息。关于此部分理论推导详见（《有限元基础教程》曾攀编著），书中详细介绍了推导过程，并且附带Matlab编程代码和Ansys实例分析。

Bar2D2Node.h函数声明文件

//***************************计算函数***************************//
/*
函数名称:		创建刚度矩阵（全局）
*/
virtual Matrix	CreateK();/*
函数名称:		计算杆单元应变
*/
virtual double	CalEpsilon();/*
函数名称:		计算杆单元应力
*/
virtual double	CalSigama();

Bar2D2Node.cpp函数实现文件

//***************************计算函数***************************//
//创建全局刚度矩阵
Matrix Bar2D2Node::CreateK()
{//计算方向正弦余弦double dealtaX = this->m_Point1->GetX() - this->m_Point0->GetX();double dealtaY = this->m_Point1->GetY() - this->m_Point0->GetY();double barLength = sqrt(pow(dealtaX ,2) + pow(dealtaY, 2));double cosTheat = dealtaX / barLength;double sinTheat = dealtaY / barLength;//构建旋转矩阵 double T[8] = {cosTheat , sinTheat, 0, 0, 0, 0, cosTheat, sinTheat};Matrix MatT(2, 4, T);//局部刚度矩阵double k[4] = {1.0, -1.0, -1.0, 1.0};Matrix Matk(2, 2, k);//刚度矩阵系数double var1 = this->m_E * this->m_A / barLength;Matk = Matk.MultNum(var1);//全局刚度矩阵this->m_MatK = MatT.Transpose().MultMat(Matk).MultMat(MatT);return this->m_MatK;
}//计算杆单元应变
double Bar2D2Node::CalEpsilon()
{//计算方向正弦余弦double dealtaX = this->m_Point1->GetX() - this->m_Point0->GetX();double dealtaY = this->m_Point1->GetY() - this->m_Point0->GetY();double barLength = sqrt(pow(dealtaX, 2) + pow(dealtaY, 2));double cosTheat = dealtaX / barLength;double sinTheat = dealtaY / barLength;//构建旋转矩阵 double T[8] = { cosTheat , sinTheat, 0, 0, 0, 0, cosTheat, sinTheat };Matrix MatT(2, 4, T);//构建B矩阵double B[2] = { -1 / barLength, 1 / barLength };Matrix MatB(1, 2, B);//构建位移矩阵double q[4] = { this->m_Point0->GetU(), this->m_Point0->GetV(), this->m_Point1->GetU(), this->m_Point1->GetV() };Matrix Matq(4, 1, q);//计算Epsilonthis->m_Epsilon = MatB.MultMat(MatT).MultMat(Matq).GetMatrixEle(0, 0);return this->m_Epsilon;
}//计算杆单元应力
double Bar2D2Node::CalSigama()
{//根据本构方程 计算sigamathis->m_Sigama = this->GetE() * this->GetEpsilon();return this->m_Sigama;
}

1.2、protected variable

保护类型变量，在后续中Bar2D2Node类会成为其他类型单元的父类，比如三维杆单元（Bar3D2Node）、二维梁单元（Beam2D2Node）等等，有些成员变量会在所派生出来的类中进行调用，此类变量都放在protected类型中。

protected:Matrix		m_MatK;					//杆单元刚度矩阵 全局double		m_Sigama;				//杆单元应力double		m_Epsilon;				//杆单元应变

1.3、private variable

私有变量，只能在类内进行调用。

private:int			ID;						//杆单元编号Point2D*	m_Point0;				//杆单元节点1号，存在顺序关系Point2D*	m_Point1;				//杆单元节点2号，存在顺序关系double		m_E;					//弹性模量double		m_A;					//杆单元横截面积

1.4、全部源码

Bar2D2Node.h函数声明文件

#ifndef _BAR_2D_H
#define _BAR_2D_H#include <iostream>
#include <list>
#include <math.h>
#include "Point.h"
#include "Matrix.h"class Bar2D2Node
{
public://***********************构造函数析构函数***********************///*函数名称:		无参构造函数*/Bar2D2Node();/*函数名称:		有参构造函数id:				ID*p0:			起始点*p1:			结束点E:				弹性模量A:				横截面积*/Bar2D2Node(int id, Point2D* p0, Point2D* p1, double E, double A);/*函数名称:		析构函数*/~Bar2D2Node();//***********************设置数值接口函数***********************///*函数名称:		设置Bar2D2Node杆单元id:				ID*p0:			起始点*p1:			结束点E:				弹性模量A:				横截面积*/virtual void		SetBar(int id, Point2D* p0, Point2D* p1, double E, double A);/*函数名称:		设置Bar2D2Node起始点与结束点*p0:			起始点*p1:			结束点*/virtual void		SetPoint(Point2D* p0, Point2D* p1);/*函数名称:		设置Bar2D2Node IDid:				ID*/void		SetID(int id);/*函数名称:		设置Bar2D2Node弹性模量E:				弹性模量*/void		SetE(double E);/*函数名称:		设置Bar2D2Node截面面积A:				横截面积*/void		SetA(double A);//***********************获取数值接口函数***********************///*函数名称:		获取Bar2D2Node起始点*/Point2D*	GetPoint0();/*函数名称:		获取Bar2D2Node终止点*/Point2D*	GetPoint1();/*函数名称:		获取Bar2D2Node单元ID*/int			GetID();/*函数名称:		获取Bar2D2Node单元弹性模量*/double		GetE();/*函数名称:		获取Bar2D2Node单元面积*/double		GetA();/*函数名称:		获取Bar2D2Node单元刚度矩阵*/Matrix		GetK();/*函数名称:		获取Bar2D2Node单元应变*/double		GetEpsilon();/*函数名称:		获取Bar2D2Node单元应力*/double		GetSigama();//***************************计算函数***************************///*函数名称:		创建刚度矩阵（全局）*/virtual Matrix	CreateK();/*函数名称:		计算杆单元应变*/virtual double	CalEpsilon();/*函数名称:		计算杆单元应力*/virtual double	CalSigama();protected:Matrix		m_MatK;					//杆单元刚度矩阵 全局double		m_Sigama;				//杆单元应力double		m_Epsilon;				//杆单元应变private:int			ID;						//杆单元编号Point2D*	m_Point0;				//杆单元节点1号，存在顺序关系Point2D*	m_Point1;				//杆单元节点2号，存在顺序关系double		m_E;					//弹性模量double		m_A;					//杆单元横截面积};#endif

Bar2D2Node.cpp函数实现文件

#include "Bar2D2Node.h"//************************************bar-2D*********************************//
//无参构造函数
Bar2D2Node::Bar2D2Node()
{
}//有参构造函数
Bar2D2Node::Bar2D2Node(int id, Point2D* p1, Point2D* p2, double E, double A)
{this->ID = id;this->m_Point0 = p1;this->m_Point1 = p2;this->m_E = E;this->m_A = A;
}//析构函数
Bar2D2Node::~Bar2D2Node()
{
}//设置Bar2D2Node杆单元
void Bar2D2Node::SetBar(int id, Point2D* p0, Point2D* p1, double E, double A)
{this->ID = id;this->m_Point0 = p0;this->m_Point1 = p1;this->m_E = E;this->m_A = A;
}//设置Bar2D2Node起始点与结束点
void Bar2D2Node::SetPoint(Point2D* p1, Point2D* p2)
{this->m_Point0 = p1;this->m_Point1 = p2;
}//设置Bar2D2Node ID
void Bar2D2Node::SetID(int id)
{this->ID = id;
}//设置Bar2D2Node弹性模量
void Bar2D2Node::SetE(double E)
{this->m_E = E;
}//设置Bar2D2Node截面面积
void Bar2D2Node::SetA(double A)
{this->m_A = A;
}//获取Bar2D2Node起始点
Point2D* Bar2D2Node::GetPoint0()
{return this->m_Point0;
}//获取Bar2D2Node终止点
Point2D* Bar2D2Node::GetPoint1()
{return this->m_Point1;
}//获取Bar2D2Node单元ID
int Bar2D2Node::GetID()
{return this->ID;
}//获取Bar2D2Node单元弹性模量
double Bar2D2Node::GetE()
{return this->m_E;
}//获取Bar2D2Node单元面积
double Bar2D2Node::GetA()
{return this->m_A;
}//获取Bar2D2Node单元刚度矩阵
Matrix Bar2D2Node::GetK()
{return this->m_MatK;
}//获取Bar2D2Node单元应变
double Bar2D2Node::GetEpsilon()
{return this->m_Epsilon;
}//获取Bar2D2Node单元应力
double Bar2D2Node::GetSigama()
{return this->m_Sigama;
}//***************************计算函数***************************//
//创建全局刚度矩阵
Matrix Bar2D2Node::CreateK()
{//计算方向正弦余弦double dealtaX = this->m_Point1->GetX() - this->m_Point0->GetX();double dealtaY = this->m_Point1->GetY() - this->m_Point0->GetY();double barLength = sqrt(pow(dealtaX ,2) + pow(dealtaY, 2));double cosTheat = dealtaX / barLength;double sinTheat = dealtaY / barLength;//构建旋转矩阵 double T[8] = {cosTheat , sinTheat, 0, 0, 0, 0, cosTheat, sinTheat};Matrix MatT(2, 4, T);//局部刚度矩阵double k[4] = {1.0, -1.0, -1.0, 1.0};Matrix Matk(2, 2, k);//刚度矩阵系数double var1 = this->m_E * this->m_A / barLength;Matk = Matk.MultNum(var1);//全局刚度矩阵this->m_MatK = MatT.Transpose().MultMat(Matk).MultMat(MatT);return this->m_MatK;
}//计算杆单元应变
double Bar2D2Node::CalEpsilon()
{//计算方向正弦余弦double dealtaX = this->m_Point1->GetX() - this->m_Point0->GetX();double dealtaY = this->m_Point1->GetY() - this->m_Point0->GetY();double barLength = sqrt(pow(dealtaX, 2) + pow(dealtaY, 2));double cosTheat = dealtaX / barLength;double sinTheat = dealtaY / barLength;//构建旋转矩阵 double T[8] = { cosTheat , sinTheat, 0, 0, 0, 0, cosTheat, sinTheat };Matrix MatT(2, 4, T);//构建B矩阵double B[2] = { -1 / barLength, 1 / barLength };Matrix MatB(1, 2, B);//构建位移矩阵double q[4] = { this->m_Point0->GetU(), this->m_Point0->GetV(), this->m_Point1->GetU(), this->m_Point1->GetV() };Matrix Matq(4, 1, q);//计算Epsilonthis->m_Epsilon = MatB.MultMat(MatT).MultMat(Matq).GetMatrixEle(0, 0);return this->m_Epsilon;
}//计算杆单元应力
double Bar2D2Node::CalSigama()
{//根据本构方程 计算sigamathis->m_Sigama = this->GetE() * this->GetEpsilon();return this->m_Sigama;
}