前言
在有限元求解完成后，需要对结果进行可视化，比如位移云图、应力云图和反力云图等。在前面一篇文章中已经实现了位移动画云图可视化，本文章将继续使用MATLAB对二维四边形单元进行应力动画云图可视化。与前面位移动画云图可视化主要区别体现在数据处理上，首先要将高斯点上的应力进行外推，再对节点应力进行平均，最后得到每个节点的应力。
数据准备
第一步还是准备数据，节点和单元与前面文章相同，不同的是不仅需要准备多个时刻的节点位移数据，还要准备多个时刻的高斯点应力数据，下面提供了5个时刻的实例数据。
%% 准备数据
% 坐标
XYZ = [ 0 0;
1 0;
2 0;
0 1;
1 1;
2 1;
0 2;
1 2;
2 2;];
% 单元
EL = [ 1 2 5 4;
2 3 6 5;
4 5 8 7;
5 6 9 8;];
% 水平方向位移时刻t=0.00
DX0 = zeros(9,1);
% 水平方向位移时刻t=0.25
DX1 = [ 0; 0; 0; 0.0174337; 0.0151149; 0.0186553; 0.0318029; 0.0311055; 0.0330584;];
% 水平方向位移时刻t=0.50
DX2 = [ 0; 0; 0; 0.0348675; 0.0302299; 0.0373107; 0.0636058; 0.0622109; 0.0661167;];
% 水平方向位移时刻t=0.75
DX3 = [ 0; 0; 0; 0.0523012; 0.0453448; 0.055966; 0.0954087; 0.0933164; 0.0991751;];
% 水平方向位移时刻t=1.00
DX4 = [ 0; 0; 0; 0.0697349; 0.0604598; 0.0746214; 0.127212; 0.124422; 0.132233;];
% 所有时刻的水平位移
DX = [DX0, DX1, DX2, DX3, DX4];
% 竖直方向位移时刻t=0.00
DY0 = zeros(9,1);
% 竖直方向位移时刻t=0.25
DY1 = [ 0; 0; 0; -0.00104918; -0.0108555; -0.0209396; -0.00926313; -0.0219317; -0.0346498;];
% 竖直方向位移时刻t=0.50
DY2 = [ 0; 0; 0; -0.00209836; -0.021711; -0.0418792; -0.0185263; -0.0438634; -0.0692995;];
% 竖直方向位移时刻t=0.75
DY3 = [ 0; 0; 0; -0.00314755; -0.0325666; -0.0628188; -0.0277894; -0.0657951; -0.103949;];
% 竖直方向位移t=1
DY4 = [ 0; 0; 0; -0.00419673; -0.0434221; -0.0837584; -0.0370525; -0.0877267; -0.138599;];
% 所有时刻的竖直位移
DY = [DY0, DY1, DY2, DY3, DY4];
% 加载应力数据
load("Sigma.mat")
其中Sigma.mat为多个时刻的高斯点应力数据，由于实例为平面问题，只存在四个应力分量，依次为
σ
11
\sigma_{11}
σ11、
σ
22
\sigma_{22}
σ22、
σ
33
\sigma_{33}
σ33和
σ
12
\sigma_{12}
σ12。1-4列表示第一个单元的四个高斯点应力，5-8列表示第二个单元的四个高斯点应力，以此类推。部分数据如下(数据的获取请查看文章顶端，或者点击下载)：
数据处理
高斯点应力外推
在有限元计算中，通过在高斯点上计算单元应力，但是在后处理中需要的是节点应力，这时就需要将高斯点上的应力外推得到节点上的应力。
外推公式如下：
(
σ
1
∗
σ
2
∗
σ
3
∗
σ
4
∗
)
=
(
a
b
c
b
b
a
b
c
c
b
a
b
b
c
b
a
)
(
σ
~
I
σ
~
I
I
σ
~
I
I
I
σ
~
I
V
)
\left( \begin{array}{c} \sigma _{1}^{*}\\ \sigma _{2}^{*}\\ \sigma _{3}^{*}\\ \sigma _{4}^{*}\\ \end{array} \right) =\left( \begin{matrix} a& b& c& b\\ b& a& b& c\\ c& b& a& b\\ b& c& b& a\\ \end{matrix} \right) \left( \begin{array}{c} \tilde{\sigma}_{\mathrm{I}}\\ \tilde{\sigma}_{\mathrm{II}}\\ \tilde{\sigma}_{\mathrm{III}}\\ \tilde{\sigma}_{\mathrm{IV}}\\ \end{array} \right)
σ1∗σ2∗σ3∗σ4∗
=
abcbbabccbabbcba
σ~Iσ~IIσ~IIIσ~IV
其中
a
=
1
+
3
2
,
b
=
−
1
2
,
c
=
1
−
3
2
a=1+\frac{\sqrt{3}}{2}, b=-\frac{1}{2}\,\,, c=1-\frac{\sqrt{3}}{2}
a=1+23,b=−21,c=1−23，
σ
~
\tilde{\sigma}
σ~表示高斯点上的应力，
σ
∗
\sigma^{*}
σ∗表示节点应力。
利用上面的这个公式即可由高斯点上的应力外推出节点上的应力，但是要求高斯点应力顺序必须是由左下角开始，逆时针绕一周，如下图所示。
绕节点应力平均
进行了高斯点应力外推后会出现一种情况就是，多个高斯点外推得到同一个节点的应力，但是这个应力不相等，这时就可以将该节点的多个应力进行平均化，即
σ
←
σ
+
σ
0
2
\sigma \gets \frac{\sigma +\sigma _0}{2}
σ←2σ+σ0
其中
σ
\sigma
σ表示当前节点的应力，
σ
0
\sigma _0
σ0表示新增加的节点应力。
节点应力计算代码
下面这个MATLAB函数就是按照上面的高斯点应力外推和绕节点应力平均写出来的，并且函数中将四个应力分量补齐成了六个应力分量，便于统一计算Mises应力。
function SigmaN = getNodeStress(Sigma, XYZ, EL)
% 由高斯点应力计算节点应力, 双线性四边形单元
% Input:
% Sigma 应力
% XYZ 节点坐标
% EL 单元节点
%
% Output:
% SigmaN 节点应力
NG = size(Sigma,2)/size(EL,1); % 高斯点个数
NS = size(Sigma,1); % 应力分量个数
NE = size(EL,1); % 单元个数
% 转化矩阵(高斯点转化到节点)
a = 1 + sqrt(3)/2;
b = -1/2;
c = 1 - sqrt(3)/2;
T = [a b c b;
b a b c;
c b a b;
b c b a;];
% 高斯节点顺序(正确顺序为左下角开始逆时针循环一圈, 若不是则需要修正顺序)
Seq = [1 2 3 4];
% 计算节点应力
SigmaN = zeros(size(XYZ,1),6)*nan; % 每个节点的应力
for i = 1:NE
SigmaE = Sigma(:,(i-1)*NG+1:i*NG)'; % 单元高斯点应力
NG0 = 4;
SigmaE = SigmaE(Seq,:); % 取出高斯点应力
SigmaE = T*SigmaE; % 由高斯点应力外插得到节点应力
% 将应力分量补齐
if NS < 6
SigmaE = [SigmaE,zeros(NG0,6-NS)];
end
% 绕节点应力平均
EL0 = EL(i,:); % 当前单元节点
Sigma0 = SigmaN(EL0,:); % 取出已经储存的节点应力
idx = ~isnan(Sigma0(:,1)); % 存在应力的节点
Sigma0(idx,:) = (Sigma0(idx,:) + SigmaE(idx,:))/2; % 节点平均均值
idx0 = ~idx; % 不存在应力的节点
Sigma0(idx0,:) = SigmaE(idx0,:); % 不存在应力则直接等于节点应力
SigmaN(EL0,:) = Sigma0; % 保存
end
end
计算Mises应力
可以通过下面的公式计算Mises应力
σ
M
=
1
2
[
(
σ
11
−
σ
22
)
2
+
(
σ
22
−
σ
33
)
2
+
(
σ
33
−
σ
11
)
2
+
6
(
σ
12
2
+
σ
23
2
+
σ
31
2
)
]
\sigma _{\mathrm{M}}=\sqrt{\frac{1}{2}\left[ (\sigma _{11}-\sigma _{22})^2+(\sigma _{22}-\sigma _{33})^2+(\sigma _{33}-\sigma _{11})^2+6\left( \sigma _{12}^{2}+\sigma _{23}^{2}+\sigma _{31}^{2} \right) \right]}
σM=21[(σ11−σ22)2+(σ22−σ33)2+(σ33−σ11)2+6(σ122+σ232+σ312)]
部分代码如下，下面代码计算了每个时刻的应力，同时还需要计算所有时刻Miese应力的最大、最小值，便于绘制颜色栏。
% 加载应力数据
load("Sigma.mat")
MisesSigma = zeros(size(DX)); % 初始化每个时刻的Mises应力
NT = size(DX,2); % 云图动画帧数
for i = 1:NT
% 计算节点应力
SigmaN = getNodeStress(Sigma{i}, XYZ, EL);
% 计算Mises应力
MisesSigma(:,i) = (1/sqrt(2))*sqrt( (SigmaN(:,1)-SigmaN(:,2)).^2 + (SigmaN(:,1)-SigmaN(:,3)).^2 ...
+ (SigmaN(:,2)-SigmaN(:,3)).^2 + 6*(SigmaN(:,4).^2 + SigmaN(:,5).^2 + SigmaN(:,6).^2) );
end
% Mises应力的最大最小值
Smin = min(min(MisesSigma));
Smax = max(max(MisesSigma));
绘图
绘制应力动画云图与绘制位移动画云图基本相同，具体细节可以查看上一篇文章。不同的是patch函数中的FaceVertexCData将位移值换成了Mises应力值，设置颜色栏范围caxis([Smin,Smax])中将最大、最小位移值换成了最大、最小Mises应力值。
%% 绘图
f = figure('Name','应力云图');
hold on
width = 550; %宽度，像素数
height = 400; %高度
left = 200; %距屏幕左下角水平距离
bottem = 200; %距屏幕左下角垂直距离
f.Position = [left,bottem,width,height];
for k = 1:NT
cla % 清除上一次的绘图
colorbar('off') % 关闭上一次的颜色栏
XYZ1 = XYZ + [DX(:,k),DY(:,k)]; % 更新坐标
patch('Faces',EL,'Vertices',XYZ1,'FaceVertexCData',MisesSigma(:,k),'FaceColor','interp','EdgeColor','k');
axis tight % 坐标轴紧凑格式
axis off % 关闭坐标轴显示
% 固定横纵坐标范围
xlim([Xmin,Xmax*1.1])
ylim([Ymin,Ymax*1.0])
% 将图向左移动, 为绘制颜色栏做准备
ax = gca;
PositionNew = get(ax,'Position');
PositionNew(1) = 0.08; % 更改图形距离左边界的距离
set(gca,'Position',PositionNew);
% 绘制颜色栏
Co = colorbar('eastoutside'); % 设置颜色栏位移图像外侧
colormap(Colormap) % 设置颜色栏
set(gca, 'FontName','Times New Roman', 'FontSize',14) % 设置字体字号
Co.Label.String = 'S, Mises'; % 设置颜色栏标题
caxis([Smin,Smax]); % 设置颜色栏的最大最小范围
% 设置颜色栏显示标签个数
NTicks = 13; % 标签个数
[cmin, cmax] = caxis;
Co.Ticks = linspace(cmin,cmax,NTicks)';
% 设置颜色栏的刻度标签为科学计数法
yt = Co.Ticks; % 获取刻度位置
yt_labels = arrayfun(@(x) sprintf('%.3e', x), yt, 'UniformOutput', false); % 将刻度位置转换为科学计数法格式的字符串
Co.YTickLabel = yt_labels; % 设置刻度标签
% 让颜色栏保持同一个位置
if k == 1
colorbarPosition = get(Co,'Position');
else
Co.Position = colorbarPosition;
end
pause(0.5) % 暂停
% 保存动态图
frame=getframe(gcf);
imind=frame2im(frame);
[imind,cm] = rgb2ind(imind,256);
DelayTime = 0.5; % 单帧播放时间
if k == 1
imwrite(imind,cm,'应力动画云图.gif','gif', 'Loopcount',inf,'DelayTime',DelayTime);
else
imwrite(imind,cm,'应力动画云图.gif','gif','WriteMode','append','DelayTime',DelayTime);
end
end
MATLAB可视化结果
Abaqus可视化结果
可以看到，MATLAB可视化结果与Abaqus可视化结果基本一致。
总代码
可视化的总代码如下，制作不易，别忘了关注和点赞喔。
clc;clear;close all
set(0,'defaultfigurecolor','w');
%% 准备数据
% 坐标
XYZ = [ 0 0;
1 0;
2 0;
0 1;
1 1;
2 1;
0 2;
1 2;
2 2;];
% 单元
EL = [ 1 2 5 4;
2 3 6 5;
4 5 8 7;
5 6 9 8;];
% 水平方向位移时刻t=0.00
DX0 = zeros(9,1);
% 水平方向位移时刻t=0.25
DX1 = [ 0; 0; 0; 0.0174337; 0.0151149; 0.0186553; 0.0318029; 0.0311055; 0.0330584;];
% 水平方向位移时刻t=0.50
DX2 = [ 0; 0; 0; 0.0348675; 0.0302299; 0.0373107; 0.0636058; 0.0622109; 0.0661167;];
% 水平方向位移时刻t=0.75
DX3 = [ 0; 0; 0; 0.0523012; 0.0453448; 0.055966; 0.0954087; 0.0933164; 0.0991751;];
% 水平方向位移时刻t=1.00
DX4 = [ 0; 0; 0; 0.0697349; 0.0604598; 0.0746214; 0.127212; 0.124422; 0.132233;];
% 所有时刻的水平位移
DX = [DX0, DX1, DX2, DX3, DX4];
% 竖直方向位移时刻t=0.00
DY0 = zeros(9,1);
% 竖直方向位移时刻t=0.25
DY1 = [ 0; 0; 0; -0.00104918; -0.0108555; -0.0209396; -0.00926313; -0.0219317; -0.0346498;];
% 竖直方向位移时刻t=0.50
DY2 = [ 0; 0; 0; -0.00209836; -0.021711; -0.0418792; -0.0185263; -0.0438634; -0.0692995;];
% 竖直方向位移时刻t=0.75
DY3 = [ 0; 0; 0; -0.00314755; -0.0325666; -0.0628188; -0.0277894; -0.0657951; -0.103949;];
% 竖直方向位移t=1
DY4 = [ 0; 0; 0; -0.00419673; -0.0434221; -0.0837584; -0.0370525; -0.0877267; -0.138599;];
% 所有时刻的竖直位移
DY = [DY0, DY1, DY2, DY3, DY4];
% 加载应力数据
load("Sigma.mat")
MisesSigma = zeros(size(DX)); % 初始化每个时刻的Mises应力
NT = size(DX,2); % 云图动画帧数
for i = 1:NT
% 计算节点应力
SigmaN = getNodeStress(Sigma{i}, XYZ, EL);
% 计算Mises应力
MisesSigma(:,i) = (1/sqrt(2))*sqrt( (SigmaN(:,1)-SigmaN(:,2)).^2 + (SigmaN(:,1)-SigmaN(:,3)).^2 ...
+ (SigmaN(:,2)-SigmaN(:,3)).^2 + 6*(SigmaN(:,4).^2 + SigmaN(:,5).^2 + SigmaN(:,6).^2) );
end
% Mises应力的最大最小值
Smin = min(min(MisesSigma));
Smax = max(max(MisesSigma));
% 位置的最大最小值
Xmin = min(min(XYZ(:,1)));
Xmax = max(max(XYZ(:,1)));
Ymin = min(min(XYZ(:,2)));
Ymax = max(max(XYZ(:,2)));
% 颜色栏
Colormap = [0 0 1;
0 0.364705882000000 1;
0 0.725490196000000 1;
0 1 0.909803922000000;
0 1 0.545098039000000;
0 1 0.180392157000000;
0.180392157000000 1 0;
0.545098039000000 1 0;
0.909803922000000 1 0;
1 0.725490196000000 0;
1 0.364705882000000 0;
1 0 0;];
%% 绘图
f = figure('Name','应力云图');
hold on
width = 550; %宽度，像素数
height = 400; %高度
left = 200; %距屏幕左下角水平距离
bottem = 200; %距屏幕左下角垂直距离
f.Position = [left,bottem,width,height];
for k = 1:NT
cla % 清除上一次的绘图
colorbar('off') % 关闭上一次的颜色栏
XYZ1 = XYZ + [DX(:,k),DY(:,k)]; % 更新坐标
patch('Faces',EL,'Vertices',XYZ1,'FaceVertexCData',MisesSigma(:,k),'FaceColor','interp','EdgeColor','k');
axis tight % 坐标轴紧凑格式
axis off % 关闭坐标轴显示
% 固定横纵坐标范围
xlim([Xmin,Xmax*1.1])
ylim([Ymin,Ymax*1.0])
% 将图向左移动, 为绘制颜色栏做准备
ax = gca;
PositionNew = get(ax,'Position');
PositionNew(1) = 0.08; % 更改图形距离左边界的距离
set(gca,'Position',PositionNew);
% 绘制颜色栏
Co = colorbar('eastoutside'); % 设置颜色栏位移图像外侧
colormap(Colormap) % 设置颜色栏
set(gca, 'FontName','Times New Roman', 'FontSize',14) % 设置字体字号
Co.Label.String = 'S, Mises'; % 设置颜色栏标题
caxis([Smin,Smax]); % 设置颜色栏的最大最小范围
% 设置颜色栏显示标签个数
NTicks = 13; % 标签个数
[cmin, cmax] = caxis;
Co.Ticks = linspace(cmin,cmax,NTicks)';
% 设置颜色栏的刻度标签为科学计数法
yt = Co.Ticks; % 获取刻度位置
yt_labels = arrayfun(@(x) sprintf('%.3e', x), yt, 'UniformOutput', false); % 将刻度位置转换为科学计数法格式的字符串
Co.YTickLabel = yt_labels; % 设置刻度标签
% 让颜色栏保持同一个位置
if k == 1
colorbarPosition = get(Co,'Position');
else
Co.Position = colorbarPosition;
end
pause(0.5) % 暂停
% 保存动态图
frame=getframe(gcf);
imind=frame2im(frame);
[imind,cm] = rgb2ind(imind,256);
DelayTime = 0.5; % 单帧播放时间
if k == 1
imwrite(imind,cm,'应力动画云图.gif','gif', 'Loopcount',inf,'DelayTime',DelayTime);
else
imwrite(imind,cm,'应力动画云图.gif','gif','WriteMode','append','DelayTime',DelayTime);
end
end
function SigmaN = getNodeStress(Sigma, XYZ, EL)
% 由高斯点应力计算节点应力, 双线性四边形单元
% Input:
% Sigma 应力
% XYZ 节点坐标
% EL 单元节点
%
% Output:
% SigmaN 节点应力
NG = size(Sigma,2)/size(EL,1); % 高斯点个数
NS = size(Sigma,1); % 应力分量个数
NE = size(EL,1); % 单元个数
% 转化矩阵(高斯点转化到节点)
a = 1 + sqrt(3)/2;
b = -1/2;
c = 1 - sqrt(3)/2;
T = [a b c b;
b a b c;
c b a b;
b c b a;];
% 高斯节点顺序(正确顺序为左下角开始逆时针循环一圈, 若不是则需要修正顺序)
Seq = [1 2 3 4];
% 计算节点应力
SigmaN = zeros(size(XYZ,1),6)*nan; % 每个节点的应力
for i = 1:NE
SigmaE = Sigma(:,(i-1)*NG+1:i*NG)'; % 单元高斯点应力
NG0 = 4;
SigmaE = SigmaE(Seq,:); % 取出高斯点应力
SigmaE = T*SigmaE; % 由高斯点应力外插得到节点应力
% 将应力分量补齐
if NS < 6
SigmaE = [SigmaE,zeros(NG0,6-NS)];
end
% 绕节点应力平均
EL0 = EL(i,:); % 当前单元节点
Sigma0 = SigmaN(EL0,:); % 取出已经储存的节点应力
idx = ~isnan(Sigma0(:,1)); % 存在应力的节点
Sigma0(idx,:) = (Sigma0(idx,:) + SigmaE(idx,:))/2; % 节点平均均值
idx0 = ~idx; % 不存在应力的节点
Sigma0(idx0,:) = SigmaE(idx0,:); % 不存在应力则直接等于节点应力
SigmaN(EL0,:) = Sigma0; % 保存
end
end
总结
这只是一个基础的示例，实际中还会有更具体的、更细致的要求，这就需要再做额外调整；另外本人也仍在学习中，这只是个人的学习笔记，可能还有一些不足之处，欢迎指正。
扫一扫
3460
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
