网站建设找什么工作室内蒙和城乡建设部网站
本系列文章致力于实现“手搓有限元，干翻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; //杆单元横截面积};#endifBar2D2Node.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;
}