一维梁有限元分析是一种常见的结构分析方法。在这种分析中,我们将一维梁划分为多个有限元,每个元素都有其自身的力学属性,例如刚度和质量。然后,我们可以使用这些元素的属性来分析整个结构的行为。
以下是一维梁有限元分析的基本步骤:
离散化:首先,我们将一维梁划分为一系列的有限元素。每个元素都有两个节点,每个节点都有一个或多个自由度(例如,旋转和位移)。
元素刚度矩阵:然后,我们对每个元素计算其刚度矩阵。对于一维梁元素,刚度矩阵通常是一个4x4的矩阵,表示了节点位移与内力之间的关系。
全局刚度矩阵:接下来,我们将所有元素的刚度矩阵组装成一个全局刚度矩阵。全局刚度矩阵描述了整个结构的力学行为。
应用边界条件:在全局刚度矩阵中,我们需要定义边界条件,例如固定、自由或者旋转等。
求解:最后,我们通过解线性方程组 F=KU (F是外部载荷向量,K是全局刚度矩阵,U是节点位移向量) 来求解节点位移。然后,我们可以基于这些位移计算元素的应力和应变。
结果可视化:通过可视化工具,我们可以查看位移、应力和应变等物理量的分布情况。
这就是一维梁有限元分析的基本过程。需要注意的是,这只是一个简化的过程,实际的有限元分析可能会涉及到更多的步骤和细节,例如非线性分析、动态分析等。
完整Calfem3.6工具箱代码见:https://download.csdn.net/download/corn1949/88817170
一维梁模型:
我们构造两根一维梁,共有3个节点,每个节点有2个自由度(y方向位移,面内旋转),构造的全局自由度矩阵就是
Edof=[1 2*1-1 2*1 2*2-1 2*2;%梁1有两个元素 1,2 这里的2*1-1, 2*1, 2*2-1等是因为每个节点有2个自由度(y方向, 面内转角),所以全局自由度编号是节点编号的两倍
2 2*2-1 2*2 2*3-1 2*3];% 梁2有两个元素 2,3
MATLAB代码:
% 简支梁的分析
clc;close all;clear all;warning off;%清除变量
rand('seed', 100);
randn('seed', 100);
format long g;
addpath(genpath('calfem3.6'));
%-----拓扑结构 -------------------------------------------------
% 一维梁的元素的自由度=y方向移动, 面内转角
Edof=[1 2*1-1 2*1 2*2-1 2*2;%梁1有两个元素 1,2 这里的2*1-1, 2*1, 2*2-1等是因为每个节点有2个自由度(y方向, 面内转角),所以全局自由度编号是节点编号的两倍
2 2*2-1 2*2 2*3-1 2*3];% 梁2有两个元素 2,3
%----- 刚度矩阵K和载荷矢量f ---------------------
K=zeros(6);
f=zeros(6,1);
f(3)=-10000;
%----- 元素刚度矩阵 ------------------------------
E=210e9;% 弹性模量
I=2510e-8;% 惯性矩
ep=[E I];% 弹性模量E和惯性矩I
ex1=[0 3];% 元素节点坐标1
ex2=[3 9];% 元素节点坐标2
% fe单元荷载矢量(4 x 1) Ke梁刚度矩阵(4 x 4)
Ke1=beam1e(ex1,ep)
Ke2=beam1e(ex2,ep)
%----- 把两根梁组装起来 将Ke组装成K---------------------------------------
K=assem(Edof(1,:),K,Ke1) % 梁1
K=assem(Edof(2,:),K,Ke2) % 梁2
%-----求解方程组并计算支撑力 -
bc=[1 0;
5 0];% 边界条件矩阵 全局节点编号的1,5的位移值被约束为0, 就是不能动
[a,r]=solveq(K,f,bc)
%----- 截面力 -------------------------------------------
Ed=extract_ed(Edof,a);
[es1,edi1]=beam1s(ex1,ep,Ed(1,:),[0],7) % 计算二维梁单元1中的截面力 梁1分为4个节点
[es2,edi2]=beam1s(ex2,ep,Ed(2,:),[0],11) % 计算二维梁单元2中的截面力 梁2分为7个节点
%----- 绘制变形梁 -----------------
figure;
hold on;
plot([0 9],[0 0]);
% plot([0,0:1:3,3:1:9,9],[0;edi1(:,1);edi2(:,1);0],'ro-');
plot([0,0:0.5:3,3:0.6:9,9],[0;edi1(:,1);edi2(:,1);0],'ro-');
axis([-1 10 -0.03 0.01]);
title('位移displacements');
%----- 绘制剪力图----------------------------------
figure;
hold on;
plot([0 9],[0 0]);
plot([0,0:0.5:3,3:0.6:9,9],[0;es1(:,1);es2(:,1);0],'ro-');
axis([-1 10 -8000 5000]);
set(gca, 'YDir','reverse');
title('剪切强度shear force');
%----- 绘制力矩图----------------------------------
figure;
hold on;
plot([0 9],[0 0]);
plot([0,0:0.5:3,3:0.6:9,9],[0;es1(:,2);es2(:,2);0],'ro-');
axis([-1 10 -5000 25000]);
set(gca, 'YDir','reverse');
title('moment');
%------------------------ end -----------------------------------
rmpath(genpath('calfem3.6'));
程序结果:
Ke1 =
2342666.66666667 3514000 -2342666.66666667 3514000
3514000 7028000 -3514000 3514000
-2342666.66666667 -3514000 2342666.66666667 -3514000
3514000 3514000 -3514000 7028000
Ke2 =
292833.333333333 878500 -292833.333333333 878500
878500 3514000 -878500 1757000
-292833.333333333 -878500 292833.333333333 -878500
878500 1757000 -878500 3514000
K =
2342666.66666667 3514000 -2342666.66666667 3514000 0 0
3514000 7028000 -3514000 3514000 0 0
-2342666.66666667 -3514000 2342666.66666667 -3514000 0 0
3514000 3514000 -3514000 7028000 0 0
0 0 0 0 0 0
0 0 0 0 0 0
K =
2342666.66666667 3514000 -2342666.66666667 3514000 0 0
3514000 7028000 -3514000 3514000 0 0
-2342666.66666667 -3514000 2635500 -2635500 -292833.333333333 878500
3514000 3514000 -2635500 10542000 -878500 1757000
0 0 -292833.333333333 -878500 292833.333333333 -878500
0 0 878500 1757000 -878500 3514000
a =
0
-0.00948586605957125
-0.022766078542971
-0.0037943464238285
0
0.007588692847657
r =
6666.66666666666
-2.35016242577935e-11
1.81898940354586e-12
-1.10591744095778e-11
3333.33333333334
-2.90475976605364e-13
es1 =
-6666.66666666667 2.44051017820694e-11
-6666.66666666667 3333.33333333336
-6666.66666666667 6666.66666666669
-6666.66666666667 10000
-6666.66666666667 13333.3333333334
-6666.66666666667 16666.6666666667
-6666.66666666667 20000
edi1 =
0
-0.00471658340184237
-0.00927506903602522
-0.013517359134889
-0.0172853559307743
-0.0204209616560214
-0.022766078542971
es2 =
3333.33333333334 20000
3333.33333333334 18000
3333.33333333334 16000
3333.33333333334 14000
3333.33333333334 12000
3333.33333333334 10000
3333.33333333334 8000
3333.33333333334 6000
3333.33333333334 4000
3333.33333333334 1999.99999999999
3333.33333333334 -4.57186372093688e-12
edi2 =
-0.022766078542971
-0.0243824701195219
-0.0247694934547524
-0.0240637450199203
-0.0224018212862835
-0.0199203187250996
-0.0167558338076267
-0.0130449630051224
-0.00892430278884464
-0.00453044963005123
-1.38777878078145e-17