【有限元学习笔记】采用两点梁单元计算挠度和转角——MATLAB程序

Infww

已于 2024-08-16 02:06:36 修改

阅读量1.9k

点赞数 4

分类专栏：有限元学习笔记 MATLAB学习笔记文章标签：学习笔记 matlab 算法

于 2023-09-15 21:40:28 首次发布

本文链接：https://blog.csdn.net/qq_53350487/article/details/132909412

版权

MATLAB学习笔记同时被 2 个专栏收录

20 篇文章

订阅专栏

有限元学习笔记

4 篇文章

订阅专栏

前言

关于有限元方法中两点梁单元的简单实例，参考书为王勖成老师的《有限单元法》，如有什么不足之处，欢迎指正。

材料参数

%% 矩形截面梁参数(mm单位制, 与Abaqus的mm单位制相同)
E = 2.1e5;  % 弹性模量
v = 0.3;    % 泊松比
b = 50;   % 宽
h = 50;   % 高
l = 1000;  % 长
I = b*h^3/12;   % 截面惯性矩

单元参数

%% 单元参数
nn = 50;    % 单元个数
le = l/nn;  % 单元长度
x = 0:le:l; % 节点坐标

刚度矩阵

单元刚度矩阵

$\boldsymbol{K}^e=\frac{EI}{l^3}\left[ \begin{array}{cc:cc} 12& 6\,l& -12& 6\,l\\ 6\,l& 4\,l^2& -6\,l& 2\,l^2\\ \hdashline -12& -6\,l& 12& -6\,l\\ 6\,l& 2\,l^2& -6\,l& 4\,l^2\\ \end{array} \right] =\left[ \begin{matrix} \boldsymbol{K}_{ii}^{e}& \boldsymbol{K}_{ij}^{e}\\ \boldsymbol{K}_{ji}^{e}& \boldsymbol{K}_{jj}^{e}\\ \end{matrix} \right]$

这里的 $l$ 指单元长度。

Ke = ((E*I)/le^3) * [12    6*le    -12    6*le;
                     6*le  4*le^2  -6*le  2*le^2;
                     -12   -6*le   12     -6*le;
                     6*le  2*le^2  -6*le  4*le^2];

将单元刚度矩阵转化为分块矩阵，方便整体刚度矩阵的叠加：

KeCell = mat2cell(Ke,[2,2],[2,2]);  % 转化为分块矩阵

整体刚度矩阵

整体刚度矩阵由单元刚度矩阵按节点编号叠加而成，例如第4个单元的左、右节点编号对应为4,5，则
$\,\, \, \, \, \,\quad\quad\quad 4 \, \, \, \quad\quad\quad\quad\quad 5 \\ \boldsymbol{K}=\begin{array}{c} \\ 4\\ 5\\ \\ \end{array}\left( \begin{matrix} \cdots& \cdots& \cdots& \cdots\\ \cdots& \boldsymbol{K}_{44}^{}+\boldsymbol{K}_{44}^{e}& \boldsymbol{K}_{45}^{}+\boldsymbol{K}_{45}^{e}& \cdots\\ \cdots& \boldsymbol{K}_{54}^{}+\boldsymbol{K}_{54}^{e}& \boldsymbol{K}_{55}^{}+\boldsymbol{K}_{55}^{e}& \cdots\\ \cdots& \cdots& \cdots& \cdots\\ \end{matrix} \right)$

%% 结构刚度矩阵
K = zeros((nn+1)*2, (nn+1)*2);  % 初始化结构刚度矩阵
K = mat2cell(K,ones(1,nn+1)*2,ones(1,nn+1)*2);  % 转化为分块矩阵
% 拼接
for k = 1:nn
    for ii = 1:2
        for jj = 1:2
            K{ii+k-1,jj+k-1} = K{ii+k-1,jj+k-1} + KeCell{ii,jj};
        end
    end
end
K = cell2mat(K);    % 转化为普通矩阵

设置载荷

单元载荷向量

$\boldsymbol{P}^e=\frac{ql}{12}\left[ \begin{matrix} 6& l& 6& -l\\ \end{matrix} \right] ^{\mathrm{T}}$

%% 单元载荷
q = 1; % 均布载荷
Pe = q*le*[1/2 le/12 1/2 -le/12]';
PeCell = mat2cell(Pe,[2,2],1);

整体载荷向量

对整体刚度矩阵叠加类似，对应节点编号相加即可

%% 结构载荷
P = zeros((nn+1)*2, 1);  % 初始化结构载荷矩阵
P = mat2cell(P,ones(1,(nn+1))*2,1);  % 转化为分块矩阵
% 拼接结构载荷
for k = 1:nn
    P{k} = P{k} + PeCell{1};
    P{k+1} = P{k+1} + PeCell{2};
end
P = cell2mat(P);    % 转化为普通矩阵

位移边界条件

简支梁

在这里插入图片描述
左、右边界条件为挠度 $w = 0$ 。

悬臂梁

在这里插入图片描述
左边界条件为挠度 $w = 0$ 以及转角 $\theta=0$ 。

这两种边界条件均为0边界条件，采用“对角元素改1法”，下面代码中的boundType是为了在两种边界条件之间转化的中间变量：

%% 边界条件(对角元素改1法)
if boundType == 1
    % 位移
    K(1,:) = 0;
    K(:,1) = 0;
    K(1,1) = 1;
    P(1) = 0;
    
    % 转角
    K(2,:) = 0;
    K(:,2) = 0;
    K(2,2) = 1;
    P(2) = 0;

elseif boundType == 2
    % 左位移
    K(1,:) = 0;
    K(:,1) = 0;
    K(1,1) = 1;
    P(1) = 0;
    
    % 右位移
    K(end-1,:) = 0;
    K(:,end-1) = 0;
    K(end-1,end-1) = 1;
    P(end-1) = 0;
end

求解方程

下面求解位移方程
$K a = P$
这里的 $K$ 为整体刚度矩阵， $a$ 为广义节点位移向量， $P$ 为整体载荷向量。
整体刚度矩阵 $K$ 具有稀疏性，非零元素如下图：
在这里插入图片描述

稀疏矩阵可以采用别的快速算法进行求解（后续文章再讨论），这里先采用MATLAB自带求解器求解：

%% 求解方程组
% a为广义位移(u, theta)
a = K\P;    % MATLAB自带解法
u = a(1:2:end);
theta = a(2:2:end);

可视化

这里的绘图细节可以阅读这篇文章【MATLAB学习笔记】绘图——去除上、右边框刻度：

function Fig(X,Y,Y1,Y2,figType,name)
    set(0,'defaultfigurecolor','w');
    alpha = 0.08;
    f = figure('Name',figType);
    n = length(X);
    box off
    plot(X,Y,'ko','MarkerSize',7,'LineWidth',1.0,'MarkerIndices',2:round(n*alpha):n);
    hold on 
    plot(X,Y1,'k>','MarkerSize',7,'LineWidth',1.0,'MarkerIndices',1:round(n*alpha):n)
    plot(X,Y2,'-k','MarkerSize',5,'LineWidth',1.3)
    xlabel('Length \it l\rm/mm','Interpreter','latex')
    switch figType
        case 'displacement'
            ylabel('Deflection \it w\rm/mm','Interpreter','tex')
        case 'stress'
            ylabel('Stress \it$\sigma_x$\rm/MPa','Interpreter','latex')
    end
    le = legend('\fontname{宋体}有限元解','\fontname{Times New Roman}Abaqus\fontname{宋体}解','\fontname{宋体}解析解','Location','best');
    le.EdgeColor = 'w';

    grid on 
    ax = gca;
    ax.GridLineStyle = '--';
    ax.GridColor = [0.5,0.5,0.5];
    ax.GridAlpha = 0.5;
    set(gca,'fontname','Times New Roman','fontsize',15)
    set(gca,'XMinorTick',true)
    set(gca,'YMinorTick',true)
    set(gca,'TickLength',[0.015,0.05])
    set(gca,'lineWidth',1)

%     ax = axes('Position',get(gca,'Position'),'XAxisLocation','top','YAxisLocation','right','Color','none');
%     set(ax,'XTick', [],'YTick', []);

    % 去除上、右边框的刻度线
    box off     % 取消边框
    ax1 = axes('Position',get(gca,'Position'),'XAxisLocation','top',...
        'YAxisLocation','right','Color','none','XColor','k','YColor','k');  % 设置坐标区
    set(ax1,'linewidth',0.9)
    set(ax1,'XTick', [],'YTick', []);   % 去掉xy轴刻度
    set(ax1,'lineWidth',1)
    hold off
    exportgraphics(f,['.\Fig\',figType,name,'.png'], 'Resolution',600)
    
end

理论挠度

$\;\;\;\;\;\;\;\;\; w=\frac{q\,x\,\left( l^3-2\,l\,x^2+x^3 \right)}{24\,E\,I}$
$\;\;\;\;\;\; w=\frac{q\,x^2\,\left( 6\,l^2-4\,l\,x+x^2 \right)}{24\,E\,I}$

结果对比

下面是程序计算结果、理论结果、Abaqus结果对比图(左简支梁，右悬臂梁)：

可以看到，程序计算出来的结果与理论结果、商业软件计算结果基本一致。

总代码

注：wAbaqusXB.mat 和 wAbaqusJZ.mat 为Abaqus计算结果的数据文件，如需要请点击链接下载Abaqus悬臂梁和简支梁挠度数据，别忘了关注和点赞喔。

clc;clear;close all
%% 矩形截面梁参数(mm单位制, 与Abaqus的mm单位制相同)
E = 2.1e5;  % 弹性模量
v = 0.3;    % 泊松比
b = 50;   % 宽
h = 50;   % 高
l = 1000;  % 长
I = b*h^3/12;   % 截面惯性矩

%% 单元参数
nn = 50;    % 单元个数
le = l/nn;  % 单元长度
x = 0:le:l; % 节点坐标



%% 设置边界条件类型
% 1.悬臂梁
boundType = 1;
% % 2.简支梁
% boundType = 2;

%% 单元刚度矩阵(课本上的公式)
Ke = ((E*I)/le^3) * [12    6*le    -12    6*le;
                     6*le  4*le^2  -6*le  2*le^2;
                     -12   -6*le   12     -6*le;
                     6*le  2*le^2  -6*le  4*le^2];
KeCell = mat2cell(Ke,[2,2],[2,2]);  % 转化为分块矩阵

%% 结构刚度矩阵
K = zeros((nn+1)*2, (nn+1)*2);  % 初始化结构刚度矩阵
K = mat2cell(K,ones(1,nn+1)*2,ones(1,nn+1)*2);  % 转化为分块矩阵
% 拼接
for k = 1:nn
    for ii = 1:2
        for jj = 1:2
            K{ii+k-1,jj+k-1} = K{ii+k-1,jj+k-1} + KeCell{ii,jj};
        end
    end
end
K = cell2mat(K);    % 转化为普通矩阵

%% 单元载荷
q = 1; % 均布载荷
Pe = q*le*[1/2 le/12 1/2 -le/12]';
PeCell = mat2cell(Pe,[2,2],1);

%% 结构载荷
P = zeros((nn+1)*2, 1);  % 初始化结构载荷矩阵
P = mat2cell(P,ones(1,(nn+1))*2,1);  % 转化为分块矩阵
% 拼接结构载荷
for k = 1:nn
    P{k} = P{k} + PeCell{1};
    P{k+1} = P{k+1} + PeCell{2};
end
P = cell2mat(P);    % 转化为普通矩阵

%% 边界条件(对角元素改1法)
if boundType == 1
    % 位移
    K(1,:) = 0;
    K(:,1) = 0;
    K(1,1) = 1;
    P(1) = 0;
    
    % 转角
    K(2,:) = 0;
    K(:,2) = 0;
    K(2,2) = 1;
    P(2) = 0;

elseif boundType == 2
    % 左位移
    K(1,:) = 0;
    K(:,1) = 0;
    K(1,1) = 1;
    P(1) = 0;
    
    % 右位移
    K(end-1,:) = 0;
    K(:,end-1) = 0;
    K(end-1,end-1) = 1;
    P(end-1) = 0;
end

%% 求解方程组
% a为广义位移(u, theta)
a = K\P;    % MATLAB自带解法
u = a(1:2:end);
theta = a(2:2:end);

%% 计算理论值
if boundType == 1
    wTrueFun = @(x)q*x.^2/(24*E*I).*(x.^2 + 6*l^2 - 4*l*x); % 悬臂梁挠度方程

elseif boundType == 2
    wTrueFun = @(x)q*x/(24*E*I).*(l^3 - 2*l*x.^2 + x.^3); % 简支梁挠度方程

end
wTrue = wTrueFun(x)';

%% 加载Abaqus数据
% load sigmaAbaqusXB.mat 
load wAbaqusXB.mat 
% load sigmaAbaqusJZ.mat 
load wAbaqusJZ.mat

%% fig
figType = 'displacement';
if boundType == 1
    wAbaqus = -wAbaqusXB;
elseif boundType == 2
    wAbaqus = -wAbaqusJZ;
end
Fig(x,u,wAbaqus,wTrue,figType,num2str(boundType))   % 可视化
FigK(K) % 混淆矩阵

function Fig(X,Y,Y1,Y2,figType,name)
    set(0,'defaultfigurecolor','w');
    alpha = 0.08;
    f = figure('Name',figType);
    n = length(X);
    box off
    plot(X,Y,'ko','MarkerSize',7,'LineWidth',1.0,'MarkerIndices',2:round(n*alpha):n);
    hold on 
    plot(X,Y1,'k>','MarkerSize',7,'LineWidth',1.0,'MarkerIndices',1:round(n*alpha):n)
    plot(X,Y2,'-k','MarkerSize',5,'LineWidth',1.3)
    xlabel('Length \it l\rm/mm','Interpreter','latex')
    switch figType
        case 'displacement'
            ylabel('Deflection \it w\rm/mm','Interpreter','tex')
        case 'stress'
            ylabel('Stress \it$\sigma_x$\rm/MPa','Interpreter','latex')
    end
    le = legend('\fontname{宋体}有限元解','\fontname{Times New Roman}Abaqus\fontname{宋体}解','\fontname{宋体}解析解','Location','best');
    le.EdgeColor = 'w';

    grid on 
    ax = gca;
    ax.GridLineStyle = '--';
    ax.GridColor = [0.5,0.5,0.5];
    ax.GridAlpha = 0.5;
    set(gca,'fontname','Times New Roman','fontsize',15)
    set(gca,'XMinorTick',true)
    set(gca,'YMinorTick',true)
    set(gca,'TickLength',[0.015,0.05])
    set(gca,'lineWidth',1)

%     ax = axes('Position',get(gca,'Position'),'XAxisLocation','top','YAxisLocation','right','Color','none');
%     set(ax,'XTick', [],'YTick', []);

    % 去除上、右边框的刻度线
    box off     % 取消边框
    ax1 = axes('Position',get(gca,'Position'),'XAxisLocation','top',...
        'YAxisLocation','right','Color','none','XColor','k','YColor','k');  % 设置坐标区
    set(ax1,'linewidth',0.9)
    set(ax1,'XTick', [],'YTick', []);   % 去掉xy轴刻度
    set(ax1,'lineWidth',1)
    hold off
    exportgraphics(f,['.\Fig\',figType,name,'.png'], 'Resolution',600)
    
end

function FigK(K)
    K(K~=0) = 1;
    K = ~K;

    figure('Name','混淆矩阵')
    clims=[0 1];  % 数值小于0的为白色，大于1的为黑色,中间的均匀分布
    imagesc(K,clims);
    colormap(gca,"gray");
    set(gca,'fontname','Times New Roman','fontsize',15)
    set(gca,'XMinorTick',true)
    set(gca,'YMinorTick',true)
    set(gca,'TickLength',[0.015,0.05])
    set(gca,'lineWidth',1)
    
    xlabel('\fontname{宋体}列')
    ylabel('\fontname{宋体}行')
    exportgraphics(gca,'.\Fig\混淆矩阵.png', 'Resolution',600)

end