计算力学——有限元编程实现

最新推荐文章于 2023-06-21 17:56:49 发布
最新推荐文章于 2023-06-21 17:56:49 发布
阅读量6.9k 收藏 70
点赞数 21
分类专栏： C++ 文章标签： c++
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_39310341/article/details/115254646
版权

计算力学——有限元编程实现

本项目采用C++编写，主要实现了平面结构三角形三节点单元、四节点四边形等参元和八节点四边形等参元以及相应的节点荷载和线性荷载处理方法，但未实现网格的自动划分算法。

此代码完成于2020年寒假期间，为了实现相应功能查阅了许多资料也借鉴了很多大佬的博客，
为找回自己的学习状态同时回顾一下计算力学的知识，写下此文

为能够更形象的表示代码实现的内容，以下述具体问题为例：

（该问题为我所学课程的大作业希望老师大大不要怪罪我）
该问题为我所学课程的大作业希望老师大大不要怪罪我

问题分析

采用C++语句面向对象的方式实现目标功能，基于PPT所示框架引入类的思想，将单元视为基础类，创建一个供不同形状单元继承的接口类 “ShapeInterface”, 在接口类中声明不同类型单元共有的功能，即：计算单元应力阵、应变阵、刚度矩阵和总刚集成方法、总荷载集成方法。
让不同类型单元继承接口类并实现相应的功能以及单元特有功能，如：计算等参元的 ”Jacobi“ 矩阵等。
通过上述的方法即可实现：由单元自行求解 ”单元刚度和荷载“ 并自动集成到 ”总体刚度和荷载“ 的目标。
从而将问题分解为：”各类型单元各自求解方法的实现“ 和 ”主程序中对单元的初始化”。

UML图

下图为接口类 “ShapeInterface”, 异常类 “ ShapeFailure”, 以及三种单元的UML类图。

Tri3 : 表示三节点三角形单元
Qua4 : 表示四节点四边形等参元
Qua8 : 表示八节点四边形等参元
（因为非计算机科班因此类图的绘制可能存在许多问题，这里贴出来的目的仅为大家更好理解整个代码的结构）

在这里插入图片描述

单元初始化

根据计算力学所学到的知识，单元刚度矩阵和单元节点荷载的求解分别需要如下参数：
在这里插入图片描述
可以看出实质上只有材料属性，节点编号以及节点坐标属于单元的固有属性，而单元的节点荷载对于不同的荷载输入会有不同的响应，因此在初始化单元时仅需要刚度矩阵 ”材料属性“, ”节点编号“, ”节点坐标“三项数据；

本程序采用文件读入的方式实现数据的写入，以三角形网格划分为例，其文件部分数据如下所示：
在这里插入图片描述

单元类的实现

以最复杂的八节点四边形等参单元为例，依次说明该类声明的各方法的实现
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述

说明

由于我删除了写报告的Typora源文件，所以只好以留档的PDF截图介绍单元类的实现，希望能够帮助大家理解每一个方法具体实现了什么。

此外，由于自己实现矩阵运算太过麻烦且已经在数值分析课程的编程作业里编写过矩阵运算的代码因此在本项目中所有的矩阵运算均采用了开源矩阵计算工具：Eigen
对于Eigen的学习和使用可以参考 Eigen: C++开源矩阵计算工具——Eigen的简单用法

最后贴上 Shape.h 文件代码和整体工程文件（写的很烂又很长以至于贴上后太长了所以换成上传的方式）

设置应该是0积分，仅供参考使用哦！！！

[文件审核中，下次再贴上来…]

Shape.h代码

/***************************************************	
* Copyright(C): 本代码所代码所有版权为gzx所有									
* @brief      : 主要实现了“三节点三角形”、“四节点四边形”和“八节点四边形”	
				单元受线性荷载和集中荷载作用下的有限元分析计算					
* @autor      : gzx																
* @version    : 1.0.0.0															
* @date       : 2020/2/20     00:27												
**************************************************/
#include <iostream>
#include <Eigen/Dense>
#pragma once
using namespace Eigen;
/**
 * @brief ShapeSpace命名空间用于编写不同形状单元
 * 该命名空间实现了异常类、接口类以及继承接口类的不同形状单元
 */
namespace ShapeSpace
{
	/**
	 * @brief ShapeFailure-异常类
	 * 实现非正常输入的报错提醒
	 */
	class ShapeFailure 
	{
	private:
		const char* Message;
	public:
		ShapeFailure(const char* msg);
		const char* GetMess();
	};
	/**
	 * @brief ShapeInterface-接口类
	 * 通过虚函数方式声明不同单元共有的方法用于形状子类的继承
	 */
	class ShapeInterface
	{
	public:
		bool CalcBMtr(double s, double t)throw(ShapeFailure);											///计算B矩阵
		virtual bool CalcDMtr()throw(ShapeFailure) = 0;													///计算D矩阵
		virtual bool CalcStif()throw(ShapeFailure) = 0;													///计算K矩阵
		virtual bool CalcLoad(MatrixXd load, MatrixXd node, MatrixXd& gF)throw(ShapeFailure) = 0;		///计算F矩阵
		virtual bool AsseStif(MatrixXd& gK)throw(ShapeFailure) = 0;										///集成gK
		virtual bool CalcCooAfter(MatrixXd& Disp)throw(ShapeFailure) = 0;								///计算单元变形后坐标
		virtual MatrixXd GetCooBefore()throw(ShapeFailure) = 0;											///返回受荷前单元坐标
		virtual MatrixXd GetCooAfter()throw(ShapeFailure) = 0;											///返回受荷后单元坐标
	};
	/**
	 * @brief Tri3-三节点三角形单元
	 * 实现三节点三角形单元的刚度矩阵、荷载矩阵等的求解
	 */
	class Tri3:public ShapeInterface
	{
	private:
		MatrixXd Ele_Mat;	///单元材料属性
		MatrixXd Nod_Num;	///单元结点编号
		MatrixXd Nod_Coo;	///单元结点坐标
		MatrixXd Nod_Co2;	///变形结点坐标
		MatrixXd Ele_Stif;	///单元刚度矩阵
		MatrixXd Ele_DMtr;	///单元弹性矩阵
		MatrixXd Ele_BMtr;	///单元应变矩阵
		double Ele_Area;	///单元面积
	public:
		const int iNum = 3;
		Tri3(MatrixXd& ele_mat, MatrixXd& nod_num, MatrixXd& nod_coo);
		bool CalcBMtr()throw(ShapeFailure);
		bool CalcDMtr()throw(ShapeFailure);
		bool CalcStif()throw(ShapeFailure);
		bool CalcLoad(MatrixXd load, MatrixXd node, MatrixXd& gF)throw(ShapeFailure);
		bool AsseStif(MatrixXd& gK)throw(ShapeFailure);
		bool CalcCooAfter(MatrixXd& Disp)throw(ShapeFailure);
		MatrixXd GetCooBefore()throw(ShapeFailure);
		MatrixXd GetCooAfter()throw(ShapeFailure);
	};
	/**
	 * @brief Qua4-四节点四边形单元
	 * 实现四节点四边形单元的刚度矩阵、荷载矩阵等的求解
	 */
	class Qua4:public ShapeInterface
	{
	private:
		MatrixXd Ele_Mat;	///单元材料属性
		MatrixXd Nod_Num;	///单元结点编号
		MatrixXd Nod_Coo;	///单元结点坐标
		MatrixXd Nod_Co2;	///变形结点坐标
		MatrixXd Ele_Stif;	///单元刚度矩阵
		MatrixXd Ele_DMtr;	///单元弹性矩阵
		MatrixXd Ele_BMtr;	///单元应变矩阵
		MatrixXd GaussW;	///Gauss加权值
		MatrixXd GaussX;	///Gauss积分点
		double Ele_Area;	///单元面积
		double detJacobi;	///J矩阵行列式
	public:
		const int iNum = 4;
		Qua4(MatrixXd& ele_mat, MatrixXd& nod_num, MatrixXd& nod_coo);
		bool CalcBMtr(double s, double t)throw(ShapeFailure);
		bool CalcDMtr()throw(ShapeFailure);
		bool CalcStif()throw(ShapeFailure);
		bool CalcLoad(MatrixXd load, MatrixXd node, MatrixXd& gF)throw(ShapeFailure);
		bool AsseStif(MatrixXd& gK)throw(ShapeFailure);
		bool CalcCooAfter(MatrixXd& Disp)throw(ShapeFailure);
		MatrixXd GetCooBefore()throw(ShapeFailure);
		MatrixXd GetCooAfter()throw(ShapeFailure);
		MatrixXd CalcJacobi(MatrixXd x, MatrixXd y, MatrixXd& Ns, MatrixXd& Nt)throw(ShapeFailure);
	};
	/**
	 * @brief Qua8-八节点四边形单元
	 * 实现八节点四边形单元的刚度矩阵、荷载矩阵等的求解
	 */
	class Qua8 :public ShapeInterface
	{
	private:
		MatrixXd Ele_Mat;	///单元材料属性
		MatrixXd Nod_Num;	///单元结点编号
		MatrixXd Nod_Coo;	///单元结点坐标
		MatrixXd Nod_Co2;	///变形结点坐标
		MatrixXd Ele_Stif;	///单元刚度矩阵
		MatrixXd Ele_DMtr;	///单元弹性矩阵
		MatrixXd Ele_BMtr;	///单元应变矩阵
		MatrixXd GaussW;	///Gauss加权值
		MatrixXd GaussX;	///Gauss积分点
		double Ele_Area;	///单元面积
		double detJacobi;	///J矩阵行列式
	public:
		const int iNum = 8;
		Qua8(MatrixXd& ele_mat, MatrixXd& nod_num, MatrixXd& nod_coo);
		bool CalcBMtr(double s, double t)throw(ShapeFailure);
		bool CalcDMtr()throw(ShapeFailure);
		bool CalcStif()throw(ShapeFailure);
		bool CalcLoad(MatrixXd load, MatrixXd node, MatrixXd& gF)throw(ShapeFailure);
		bool AsseStif(MatrixXd& gK)throw(ShapeFailure);
		bool CalcCooAfter(MatrixXd& Disp)throw(ShapeFailure);
		MatrixXd GetCooBefore()throw(ShapeFailure);
		MatrixXd GetCooAfter()throw(ShapeFailure);
		MatrixXd CalcJacobi(MatrixXd x, MatrixXd y, MatrixXd& Ns, MatrixXd& Nt)throw(ShapeFailure);
	};
}

应该没人看到这了吧，写这么多好像暴露学校了，(lll￢ω￢)，不过这门课应该没什么人用C++写吧
sUaN_SuAn____
关注
  • 21
    点赞
  • 踩
  • 70
    收藏
    觉得还不错? 一键收藏
  • 5
    评论
专栏目录
有限元之编程示例三：有限元求抛物型偏微分方程
欢迎各位大佬光临“寒舍”！请多多指点，共同进步！
06-05 270
主要对有限元法求解抛物型偏微分方程边值问题进行C++编程。
用C++写一个简单的有限元程序
weixin_42584586的博客
02-13 548
下面是一个用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 1238
平面三角形单元有限元实现 C++
用c++语言完成有限元计算
weixin_42601702的博客
02-09 478
用C语言完成有限元计算需要使用到数学知识，特别是微分方程。首先，需要编写C语言程序来解决所给出的有限元问题，并计算出所需的结果。其次，需要使用有限元分析和数值分析等方法来求解有限元问题，以获得准确的解。 ...
等参元：平面四节点四边形等参元的刚度矩阵的计算
qq_39642029的博客
05-09 5613
补充： 其中： -------------------------------------------------------------------------------------- 例题： 如图为一个平面3节点四边形等参元，采用4点高斯积分计算该单元刚度矩阵。 [1有限元基础教程-曾攀 ...
稀疏矩阵的快速转制c语言代码_有限元方法：从代码入手(4)
weixin_29743937的博客
12-18 427
稀疏矩阵的封装由有限元方法获得的单元刚度矩阵通常与单个单元相关联。这样单元矩阵与有限个单元节点相互关联，由这些单元矩阵组装成的全局矩阵就含有大量的非零元，即稀疏矩阵。为了节省存储内存，通常仅仅存储非零元的信息。稀疏矩阵的具有多种存储格式，包括行主元存储，列主元存储，坐标存储格式。在这里封装时，仅仅对CSR格式的稀疏阵进行封装。以上图所示的稀疏矩阵为例，在存储的时候，需要一个名为Value的数组存储...
C++经典程序代码大全
热门推荐
zhusg的专栏
04-14 16万+
//根据半径计算圆的周长和面积 #include const float PI=3.1416; //声明常量(只读变量)PI为3.1416 float fCir_L(float); //声明自定义函数fCir_L()的原型 float fCir_S(float); //声明自定义函数fCir_S()的原型 //以下是main()函数
有限元编程作业——c语言三角形常应变单元计算程序
最新发布
qwersdf_we的博客
06-21 923
节点坐标：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。其中弹性模量设置为1，泊松比设置为0.2。边界条件：节点3,4的位移为0。
平面四边形八节点-等参元程序.zip_finite element_八节点_平面八节点_等参元
09-23
有限元的8节点四边形等参单元，带有详细说明，MATLAB平台编制。
有限元编程：基础篇
slw213106的博客
12-18 3148
、前言 相信很多做过有限差分之后又想做做有限元的初学者会有和我一样的困惑，能看懂有限元算法的理论分析，但是真正应用到实际编程当中之前心里发怵，废话不多说，求人不如求己，看懂这篇文章将会让你迅速掌握有限元最基础的编程思想。 二、以经典扩散方程为例(反常扩散方程可类比此例) 考虑如下扩散方程初边值问题 ⎧⎩⎨⎪⎪⎪⎪∂u(x,t)∂t=∂2u(x,t)∂x2+f(x,t),u(a,t)=0,u(b,t)=0,u(x,0)=φ(x),a<x<b,0<t≤T,0≤t≤T,a≤x≤b
四边形单元刚度矩阵程序.zip_单元刚度矩阵_四边形 有限元_四边形单元_四边形单元刚度矩阵
07-14
实现有限元中四边形单元的刚度矩阵计算，实现参数转换功能
计算结构力学
12-10
秦荣老师编写，很好的力学书
力学研究所计算力学
09-24
计算力学最小势能原理表述为：对于线弹性结构，满足连续性要求和位移边界条件的位移中，满足平衡方程和力边界条件的位移使势能达到最小值
有限元法及其程序实现.pdf
02-12
有限元方法的概论，翻译自国外书。附带fortran程序程序实现
C++编程有限元公式
05-13
适合于C++有限元编程初学者，利用C++编辑有限元力学公式，作用初学的参考，可以借鉴。
MATLAB有限元分析与应用——包含完整的源代码
12-27
在有限元方法的世界里，MATLAB是一个强大的工具，其易用性和灵活性使得它成为科研和工程领域进行数值计算的首选平台。本书的目标是帮助读者理解有限元分析的基本概念，并学会如何在MATLAB环境下实施这一方法。 有限...
Pinball HDU - 6373（物理受力分析？）
Sirius_han的博客
08-09 536
Pinball  题目链接：HDU - 6373 题意：一个小球有高空落到斜坡，忽略能量损失，小球能在斜坡上弹几次； 思路：高中最简单的受力分析都不会了，，，，，，； 终于将受力分析分析对了后就是上图；先算出小球落到斜面的速度，再求垂直斜面的运动周期，沿斜面方向的运动时间，无情的WA了；下面在上两张图：                 这两张图的受力分析由小球开始...
c语言编写简支梁有限元程序
06-06
简支梁有限元程序的核心是构建刚度矩阵和载荷向量，然后通过求解线性方程组得到位移和反力。以下是一个简单的C语言程序示例，可以用于求解简支梁的位移和反力： ```c #include <stdio.h> #include <stdlib.h> #include <math.h> #define N 3 //每个单元节点数 #define M 2 //梁的数量 #define DOF 2 //每个节点的自由度数 void assemble(double *K, double *f); void solve(double *K, double *f, double *u); int main() { //定义节点坐标 double x[N * M] = {0.0, 0.0, 2.0, 0.0, 4.0, 0.0, 6.0, 0.0, 8.0, 0.0, 10.0, 0.0}; //定义单元节点编号 int e[N * (M - 1)] = {1, 2, 3, 4, 5, 6}; //定义杨氏模量和截面面积 double E = 2.1e11, A = 0.01; //定义节点载荷 double f[N * M * DOF] = {0.0, 0.0, 0.0, 0.0, 0.0, -1000.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1000.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1000.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1000.0}; double K[N * M * DOF][N * M * DOF] = {0.0}; double u[N * M * DOF] = {0.0}; assemble(&K[0][0], &f[0]); solve(&K[0][0], &f[0], &u[0]); for (int i = 0; i < N * M * DOF; i++) { printf("u[%d] = %.4f\n", i, u[i]); } return 0; } void assemble(double *K, double *f) { //定义局部刚度矩阵 double k[N * DOF][N * DOF] = {0.0}; double L, c, s; for (int i = 0; i < M; i++) { //计算单元长度 L = sqrt(pow(x[N * i + 2] - x[N * i], 2.0) + pow(x[N * i + 3] - x[N * i + 1], 2.0)); //计算单元方向余弦 c = (x[N * i + 2] - x[N * i]) / L; s = (x[N * i + 3] - x[N * i + 1]) / L; //计算局部刚度矩阵 k[0][0] = k[1][1] = E * A / L; k[0][1] = k[1][0] = -E * A / L; k[2][2] = k[3][3] = 12 * E * A / pow(L, 3.0); k[2][3] = k[3][2] = 6 * E * A / pow(L, 2.0); k[0][2] = k[2][0] = k[1][3] = k[3][1] = -E * A / L; k[1][2] = k[2][1] = k[0][3] = k[3][0] = -E * A / L; k[4][4] = k[5][5] = 4 * E * A / L; k[4][5] = k[5][4] = -6 * E * A / pow(L, 2.0); k[4][0] = k[0][4] = k[5][1] = k[1][5] = -6 * E * A / pow(L, 2.0) * s; k[4][1] = k[1][4] = k[5][0] = k[0][5] = 6 * E * A / pow(L, 2.0) * c; //组装局部刚度矩阵到全局刚度矩阵 for (int j = 0; j < N; j++) { for (int k = 0; k < N; k++) { for (int m = 0; m < DOF; m++) { for (int n = 0; n < DOF; n++) { int row = j * DOF + m; int col = k * DOF + n; K[row][col] += k[m + j * DOF][n + k * DOF]; } } } } } //组装载荷向量到全局载荷向量 for (int i = 0; i < N * M * DOF; i++) { f[i] = f[i]; } } void solve(double *K, double *f, double *u) { //求解位移和反力 //对于简支梁，两个节点的横向位移是已知的，因此只需要求解纵向位移和反力 double K11 = K[2][2], K12 = K[2][5], K21 = K[5][2], K22 = K[5][5]; double f1 = f[2], f2 = f[5]; double det = K11 * K22 - K12 * K21; u[2] = (K22 * f1 - K12 * f2) / det; u[5] = (K11 * f2 - K21 * f1) / det; } ``` 该程序中，`x`数组表示节点坐标，`e`数组表示单元节点编号，`E`和`A`分别表示杨氏模量和截面面积，`f`数组表示节点载荷，`K`数组表示刚度矩阵，`u`数组表示位移和反力。程序的核心是`assemble`函数和`solve`函数，`assemble`函数用于构建刚度矩阵和载荷向量，`solve`函数用于求解线性方程组。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 5
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

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

抵扣说明：

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

余额充值