可视化一个弹性力学的解析解
过冷水诚挚邀请你加入Matlab仿真秀官方交流群进行Matlab学习、问题咨询、 Matlab相关资料下载，qq：927550334
弹性力学问题的解析解一般都是非常冗长复杂的公式，对于绝大多数学生来说，在脑海中形象化展示这些公式的结果并快速理解其含义都是不可能完成的任务。 而如果将这些公式可视化后变成直观的图形或图像， 则其含义就很容易理解了。因此，对复杂公式的可视化将非常有助于对弹性力学解析解的理解。本章以一个经典弹性力学问题——对径受压圆盘为例，来说明复杂解析解可视化的过程和方法。
9.1 对径受压圆盘的应力分布
如图9-1所示的对径受压圆盘，圆盘半径为r，厚度为t。按弹性力学理论，圆盘表面的应力解为 ：
图 9-1 巴西圆盘试件尺寸图
9.2 应力分布的可视化过程
本节以应力分量σx的可视化为例进行说明。可视化时首先要生成一个数据矩阵，这需要先在给定区域内生成一个表达σx的矩阵，然后就可以利用MATLAB的强大的绘图功能绘制各种图形。 设圆盘的半径r =25 mm，厚度为b=1 mm。绘制当p＝-1N时的应力解。
9.2.1 数据矩阵的生成
首先设定计算区域，并给出圆盘上每一点的坐标。这一步可通过MATLAB的meshgrid函数实现，其功能在上一篇已经叙述过。 已知圆盘的半径为25 mm，设定绘图区域如图9-2所示。
图 9-2 绘图区域
例1：生成如图9-2所示的绘图区域，并根据式
（9-1）求解出对应区域的σx。
>>r = 25;t = 5;p = -1; [x,y]= meshgrid(-25:0.1:25,-25:0.1:25);%%生成计算区域点阵 yrx = (y r).*(y r) x.*x; yrx = yrx.*yrx; ryx = (r-y).*(r-y) x.*x; ryx = ryx.*ryx; yr = (y r); ry = (r-y); sigmax = 2 * p / (pi*t) * (yr.*x.*x./yrx ry.*x.*x./ryx-1/(2*r));
9.2.2 绘图中的细节考虑
用上一节的方法计算出的圆盘的应力解存在一个问题，即图9-2中的圆盘之外的区域也存在应力值。为了消除这些区域的影响，可先判定圆盘的范围，然后将圆盘之外区域的数据点设为0。用一个语句可完成上述设置：
>> sigmax(x.*x y.*y - r.*r > 0)=0;
在应力分量数据矩阵生成后，就可以利用MATLAB绘图函数进行可视化了。图9-3是用contourf显示的σx的结果。可以发现尽管将圆盘之外的区域设置为0，但这些设置却给绘图函数带来“麻烦” 。在试件边界处，由于应力突跳，正式绘图区的数据很难显示出来（图9-3） 。
图 9-3 集中载荷下圆盘 σx应力分布图
（a）20 条等值线； （b）150 条等值线
利用MATLAB的nan，可以完美解决上述问题。nan有两个特点，即参与任何运算时其结果都是nan；绘图函数不处理nan数据点，在该数据点处显示背景色。
将圆盘区域之外的点设置为nan：
>> sigmax(x.*x y.*y - r.*r > 0)=nan;
将其可视化后， 可以明显地看出有用的绘图区域， 但由于端部应力梯度较大，等值线非常集中，效果不够美观，如图9-4（a）所示。可以利用shading flat命令去掉黑色的等值线，使图像更加好看些，如图9-4（b）所示。
图 9-4 巴西圆盘应力结果（σx）可视化
（a）未去掉等值线； （b）去掉等值线
9.2.3 可视化结果
也可以用MATLAB的其他绘图方式对应力解可视化显示。 例2：使用pcolor函数对巴西圆盘的x方向应力、y方向应力和剪应力进行可视化显示。
%% pcolor
figure;
h1 = subplot(1,3,1);pcolor(x,y,sigmax);caxis([-10e-3,0.001]);
axis equal;axis tight;colorbar;shading flat;
h2 = subplot(1,3,2);pcolor(x,y,sigmay);caxis([-0.1,-0.001]);
axis equal;axis tight;colorbar;shading flat;
h3 = subplot(1,3,3);pcolor(x,y,sigmaxy);caxis([-0.06,0.06]);
axis equal;axis tight;colorbar;shading flat;
title(h1,'\it\sigma_x','fontname','Times new roman');
title(h2,'\it\sigma_y','fontname','Times new roman')
title(h3,'\it\tau_{xy}','fontname','Times new roman') 采用以上方法可方便得到如图9-5所示的可视化结果。
图 9-5 巴西圆盘 σx、σy 及 τxy分布
另外，还可以通过σx、σy 及τxy计算巴西圆盘的最大剪应力τmax和主方向θ，如图9-6所示。在实验力学领域，有一种测量应力场分布的方法称为光弹方法，其等差线和等倾线正好代表的是主应力差和主方向，其条纹分布与这两个场的分布是类似的，本书后面还会详述。
tau_max = (sigmax-sigmay)/2;%计算最大剪应力
theta = sigmax; %初始化主方向
m = 1;
for i = 1:501
for j = 1:501
if ~isnan(sigmax(i,j))
[k l] = eig([sigmax(i,j),sigmaxy(i,j);sigmaxy(i,j),
sigmay(i,j)]);
if k(1,m)>=0
theta(i,j) = atan(k(2,m)/k(1,m));
elseif k(2,m)>=0
theta(i,j) = pi atan(k(2,m)/k(1,m));
elseif k(2,m)<0
theta(i,j) = pi - atan(k(2,m)/k(1,m));
end
end
end
end
figure;
h1 = subplot(1,2,1);pcolor(x,y,tau_max);caxis([0,0.02]);axis equal; ...
colorbar;shading flat;ylim([-25 25]);
h2 = subplot(1,2,2);pcolor(x,y,theta);axis equal;...
colorbar;shading flat;ylim([-25 25]);
title(h1,'{\it\tau}_{max}','fontname','Times new roman')
title(h2,'\it\theta','fontname','Times new roman')
图 9-6 巴西圆盘 τmax及主方向 θ 的分布
9.3 动画显示加载过程的应力变化
有了上述可视化程序，输入不同载荷，即可得到不同的应力场。如果想要连续显示不同载荷下应力场的变化，可以利用MATLAB的动画制作功能，生成动画来动态显示结果。
9.3.1 动画制作的基本概念
动画一般由视频文件来保存，如.avi文件。一个视频文件（不包括声轨）的基本结构如图9-7所示，除文件头外，其基本组成部分是一帧帧图像，每一帧称为frame。所以，制作动画的基本过程就是两步：一是创建视频文件，二是构建frame的内容。
图 9-7 视频文件结构
在MATLAB中创建视频文件非常简单，用函数Avifile可以方便实现，其格式为
句法
aviobj = avifile(filename, 'Param1', Val1, 'Param2', Val2,...)
说明
avifile：创建一个 AVI 格式的视频文件。
filename：视频文件名。
'Param1', Val1, 'Param2', Val2, ...：设定参数，具体可设定的参数可参照“Help” 。
创建视频文件后，就要向其中填充frame，填充完毕后，用close函数将其关闭（相当于对文件进行“封口” ）后，就可以用播放软件进行播放了。
由前面的叙述可知，生成视频文件最重要的一步是填充frame。填充frame用 addframe函数实现，其格式为：
句法
aviobj = addframe(aviobj,F)
说明
将帧图像放入视频文件中，其中 aviobj 为使用 avifile 命令创建的视频文件的对象，F为一个 frame。
一般情况下，用一个循环语句即可完成所有frame的填充。
动画制作中最有技术含量的步骤是生成frame。MATLAB中提供了两种思路
生成frame，即绘图快照式和图像转换式，下面分别介绍。
9.3.2 两种生成 frame 的方式
（1） 绘图快照式
绘图快照式是对MATLAB显示的图形或图像进行“快照” ，生成frame。其方法为：
句法
F = getframe(h)
说明
其中 h 可以为坐标轴的句柄（如“gca” ） ，也可以是窗口的句柄（如“gcf” ） 。
例3：消失的圆。用动画显示一个直径100的圆（图9-8）逐渐消失的过程。
图 9-8 消失的圆
%% 制作视频
close all;clear all;clc
mov = avifile('round_Video.avi','fps',10);
% 创建视频文件 Video.avi,设帧率为 10fps
h = figure('color',[1 1 1]);
round = rectangle('Position',[-50,-50,100,100],'Curvature',[1,1],...
'EdgeColor','k','LineWidth',4);%圆初始化
axis off;axis equal;axis tight;hold on
for i = 100:-1:0 % 设定循环
F = getframe(h);% 采集当前坐标轴内的图像
mov = addframe(mov,F); % 存储图像至视频文件
if i
set(round,'Position',[-i/2,-i/2,i,i])
else
set(round, 'FaceColor',[1 1 1],'EdgeColor',[1 1 1])
end
end % 循环结束
mov = close(mov); % 关闭视频文件（2） 图像转换式
图像转换式是专门针对图像动画而提供的一种frame生成方式，其生成用到的函数为im2frame，具体说明如下：
句法
F = im2frame(I)
说明
I：为图像的灰度矩阵。
下面以制作英文字母动态显示的视频动画为例，如图9.9所示。
图 9-9 原始图像
图9-9中包含26个灰度值各不相同的英文字母，在本实例中，需要对图像进行二值化处理，实现字母的连续显示。具体的实现程序如下。
>> close all
clear all
clc
%% 读取图像
image = imread('字母.bmp'); % 读取图像
%% 显示字母
for i = 1:26 % 设定循环,26 个字母连续显示 26 次
threshold = 1-(i*8-7)/256; % 设定灰度阈值
bw = im2bw(image,threshold); % 根据阈值对原始图像进行二值化处理
imshow(bw) % 显示图像二值化图像
end % 循环结束程序在执行过程中，会不断显示经过二值化处理后的图像，共26幅，如图9-10所示。
图 9-10 26 幅经过二值化处理后的图像
实现连续显示后就可以制作视频了。根据上文的流程介绍，只需要在源程序上添加帧图像采集和视频制作的语句即可，具体实现程序如下（加重字体的语句为添加语句） 。
>>> close all
clear all
clc
%% 读取图像
image = imread('字母.bmp'); % 读取图像
%% 制作视频
mov = avifile('Video.avi','fps',2); % 创建视频文件 Video.avi,设帧率为 2fps
for i = 1:26 % 设定循环,26 个字母连续显示 26 次
threshold = 1-(i*8-7)/256; % 设定灰度阈值
bw = im2bw(image,threshold); % 根据阈值对原始图像进行二值化处理
F = im2frame(bw); % 直接将图像转化为 frame
mov = addframe(mov,F); % 存储图像至视频文件
end % 循环结束
mov = close(mov); % 关闭视频文件9.3.3 对径受压圆盘加载过程的动画显示
利用上面的知识，可以对圆盘加载过程中x方向应力变化实现动画显示。给定圆盘的加载过程为从0.05 N至1 N，载荷步长为0.05 N。
在弹性范围内，对圆盘进行连续加载，应力与载荷成正比，因而只需将初始载荷下计算得到的应力场乘以相应的载荷倍数，即可实现连续加载。
具体程序如下：
>>> mov = avifile('bra_Video.avi','fps',10);
% 创建视频文件 Video.avi,设帧率为 10fps
h = figure('color',[1 1 1]);
for i = 0.05:0.05:1 % 设定循环
contourf(x,y,sigmax*i,500);caxis([-10e-3,0.003]);...
axis equal;axis tight;colorbar;shading flat;
F = getframe(h);% 采集当前坐标轴内的图像
mov = addframe(mov,F); % 存储图像至视频文件
end % 循环结束
mov = close(mov); % 关闭视频文件结果如图9-11所示。
图 9-11 各载荷下 σx的分布
（a）载荷 0.1 N； （b）载荷 0.5 N； （c）载荷 1.0 N
过冷水发表于仿真秀 平台原创文章，未经授权禁止私自转载，如需转载请需要和作者沟通表明授权声明，未授权文章皆视为侵权行为，必将追责。如果您希望加入Matlab仿真秀官方交流群进行Matlab学习、问题咨询、 Matlab相关资料下载均可加群：927550334。
精品回顾
过冷水和你分享 matlab读取存储各种文件的方法 文末有独家金曲分享
