挨拶 & 雑談
2022年, あけましておめでとうございます (記事を書き始めた時は正月休みでした). 本年もMATLAB Student Ambassadorをよろしくお願いいたします. さて, Twitterにおける, 2022年最初の投稿はこんな感じでした.
Happy New Year!!
— ヒビヤ@MATLAB Ambassador (@hibs_MATLAB_Amb) December 31, 2021
新年, 明けましておめでとうございます.
本年もよろしくお願いいたします.
2022年の干支にちなんで, 年末に行った有限要素法によるトラ柄の研究の成果です.#MATLAB pic.twitter.com/d18CJMWBlq
これは私のスウェーデン留学時代からの超個人的な習慣です. スウェーデンに留学していたときに構造力学の授業があり, 友達と「授業で学んだことを生かして何か作りたいね～」なんて話をしてました. その結果, こんなものが完成しました. 当時は「日本では正月になると新春とか言われるんだよー」と話しながらcolormapをspringにしてましたが, MATLAB Student Ambassadorになった2022年には, もう少し頑張ろうという事で, 2022年の干支であるトラを頑張ってカラーマップで表現してみました!!
この記事では, 使用したカラーマップに加えて, 数字を模した穴のある有限要素法解析モデルをMATALBで作るコードを紹介したいと思います!!
利用したバージョン及び製品
本記事に掲載されているコードはR2021bにおいて検証を行っています.
また, 本記事に掲載されているコードをMATLABで実行するためには, Partial Differential Equation Toolboxが必要になります.
MATLAB における colormap の作成
MATLABでは自分でcolormapを作成することができます. 方法は簡単で, 0.0~1.0の範囲の値から成る3列の行列を作成し, 1列目に赤, 2列目に緑, 3列目に緑の強度を入力します (MATLAB Documentationより). 最大値が255ではないところは要注意です!! 今回は, トラを表現するために, 以下のcolormapを作成しました.
\begin{align}
red(x) &= \left\{
\begin{array}{ll}
0 & (x \lt 0.2) \\
(5x-1)^4 & (0.2 \leq x < 0.4) \\
1 & (0.4 \leq x < 0.6) \\
(4-5x)^4 & (0.6 \leq x < 0.8) \\
0 & (x > 0.8)
\end{array}\right. \\
green(x) &= \left\{
\begin{array}{ll}
0 & (x \lt 0.2) \\
0.8431(5x-1)^4 & (0.2 \leq x < 0.4) \\
0.8431 & (0.4 \leq x < 0.6) \\
0.8431(4-5x)^4 & (0.6 \leq x < 0.8) \\
0 & (x > 0.8)
\end{array}\right. \\
blue(x) &= 0
\end{align}
n = 1000;
cmap = zeros(n,3);
color1 = [0,0,0];
color2 = [255/255,215/255,0/255];
grad_start1 = 0.2;
grad_finish1 = 0.4;
grad_start2 = 0.6;
grad_finish2 = 0.8;
order = 4;
for ii=1:n
if (ii/n < grad_start1)
cmap(ii,:) = color1;
elseif (ii/n < grad_finish1)
cmap(ii,:) = color2*(1-((grad_finish1*n-ii)/n/(grad_finish1-grad_start1))^order);
elseif (ii/n < grad_start2)
cmap(ii,:) = color2;
elseif (ii/n < grad_finish2)
cmap(ii,:) = color2*((grad_finish2*n-ii)/n/(grad_finish2-grad_start2))^order;
else
cmap(ii,:) = color1;
end
end
cmap(cmap > 1) = 1;
穴あき有限要素法モデルの作成
次に, 数字を模した穴をあけた解析モデルを作成しました. これには、本記事の最後に載せた自作関数, crate_NumModel, を利用しています. この関数では, 穴に表示する1つの整数Numと穴の大きさを指定するパラメータradius_rateを受け取ってから, 2-D Geometry Creation at Command Lineのルールに基づいた配列を作成し, decsg関数を使ってモデルを作成しています.
今回は, 2022年を祝うための解析モデルが必要なため, 以下を作成しました.
今のところ桁数に上限は設けいていないため, より大きな数字も出力が可能です.
Partial Differential Equation Toolboxを利用した有限要素法解析
この部分については, 別の記事([Blender+MATLAB] 有限要素法解析を行い, グラデーションのついた文字を作る!!)に書かれておりますので, ぜひご覧ください!!
基本的に解析において必要なのは,
- 材料定数(ヤング率, ポアソン比)の決定 (structuralProperties)
- 負荷する荷重の大きさと方向の決定 (structuralBoundaryLoad)
- 固定する部分の決定 (structuralBC)
- メッシュの作成 (generateMesh)
となります. 今回は以下のように設定の上で解析を行いました.
% Material
E = 197*10^9; % ヤング率 [Pa]
nu = 0.3; % ポアソン比
structuralProperties(model,'YoungsModulus',E,'PoissonsRatio',nu);
% Boundary Condition
structuralBoundaryLoad(model,'Edge',3,'SurfaceTraction',[100;0]);
structuralBoundaryLoad(model,'Edge',4,'SurfaceTraction',[0;100]);
structuralBC(model,'Edge',1,'XDisplacement',0);
structuralBC(model,'Edge',2,'YDisplacement',0);
% Generate Mesh
mesh = generateMesh(model,'Hmax',1/70,'Hmin',1/70,'Hgrad',2.0,'GeometricOrder','linear');
% Analysis
R = solve(model);
結果
このように, カラーマップ及び解析モデルを作成して有限要素法解析を行った結果, 見事, トラのような柄(???)の解析結果を得られました!!
figure
pdeplot(model,'XYData',R.Stress.xx,'ColorMap',cmap,'ColorBar','off')
c = colorbar('southoutside');
c.Label.String = 'Stress [MPa]';
axis equal
title 'Normal Stress Along x-Direction';
set(gca,'FontSize',14);
figure
pdeplot(model,'XYData',R.Strain.yy,'ColorMap',cmap,'ColorBar','off')
colorbar('southoutside')
axis equal
title 'Normal Strain Along y-Direction';
set(gca,'FontSize',14);
穴あきの解析モデルを作成する自作関数
R2021bでのみ検証を行っております.
function model = create_NumModel(Num,radius_rate)
%
% create_NumModel - Create FEM Model including Number
%
% model = create_NumModel(num,Radius_Rate)
%
% -Inputs-
% Num - An Integer including in the model
% Radius_Rate - Define the Size of Circle (Default : 1)
%
% -Outputs-
% model - FEM Model
%
% --- Usage Example ---
%
% model = create_NumModel(2022,1);
%
% ---------------------
%
%
switch nargin
case 1
radius_rate = 1;
end
if (Num ~= round(Num))
% Judge the Num is Integer or not
error('Error!! Number should be an Integer.');
elseif (numel(Num) > 1)
% Judge the Num is Integer or not
error('Error!! Number should be One Integer.');
else
% Sort the Input Number
numchar = char(num2str(round(Num)));
n = length(numchar);
num = zeros(1,n);
for ii =1:n
if (numchar(ii) == '0')
num(ii) = 10;
else
num(ii) = str2double(numchar(ii));
end
end
% Data of relationship between Number and Whole Location
ND = zeros(10,13);
ND( 1,:) = [0 0 1 0 1 0 0 1 0 1 0 0 1];
ND( 2,:) = [1 1 1 0 1 1 1 1 1 0 1 1 1];
ND( 3,:) = [1 1 1 0 1 1 1 1 0 1 1 1 1];
ND( 4,:) = [1 0 1 1 1 1 1 1 0 1 0 0 1];
ND( 5,:) = [1 1 1 1 0 1 1 1 0 1 1 1 1];
ND( 6,:) = [1 1 1 1 0 1 1 1 1 1 1 1 1];
ND( 7,:) = [1 1 1 1 1 0 0 1 0 1 0 0 1];
ND( 8,:) = [1 1 1 1 1 1 1 1 1 1 1 1 1];
ND( 9,:) = [1 1 1 1 1 1 1 1 0 1 1 1 1];
ND(10,:) = [1 1 1 1 1 1 0 1 1 1 1 1 1];
% Size Data
l = 1; % Plane Size
r = radius_rate*(l/4)/10; % Radius
w = 15; w2 = w/2; % Define Margin Rate
% Hole Location
hole_data = zeros(13,2);
hole_data(1,:) = [-l/w l/w2];
hole_data(2,:) = [0 l/w2];
hole_data(3,:) = [l/w l/w2];
hole_data(4,:) = [-l/w l/w];
hole_data(5,:) = [l/w l/w];
hole_data(6,:) = [-l/w 0];
hole_data(7,:) = [0 0];
hole_data(8,:) = [l/w,0];
hole_data(9,:) = [-l/w -l/w];
hole_data(10,:) = [l/w -l/w];
hole_data(11,:) = [-l/w -l/w2];
hole_data(12,:) = [0 -l/w2];
hole_data(13,:) = [l/w -l/w2];
% Create Model
model = createpde('structural','static-planestress');
R11 = [3,4,0,l*n/4,l*n/4,0,l/4,l/4,-l/4,-l/4]';
gm = R11;
sf = 'R11';
ns = "R11";
holenum = 0;
for ii=1:n
for jj=1:13
if (ND(num(ii),jj) == 1)
C1 = zeros(length(R11),1);
C1(1:4,1) = [1,(2*ii-1)*l/8+hole_data(jj,1),hole_data(jj,2),r]';
gm = [gm C1];
holenum = holenum+1;
str = sprintf('C%d',holenum);
sf = sprintf('%s-%s',sf,str);
str = sprintf("C%d",holenum);
ns = [ns;str];
end
end
end
%sf = 'R1-C1-C2-C3-C4-C5-C6-C7-C8-C9-C10-C11';
ns = char(ns);
ns = ns';
g = decsg(gm,sf,ns);
geometryFromEdges(model,g);
end
最後に
なぜ, 年始を過ぎてまで一生懸命書いたかというと, 来年の干支である卯(ウサギ)は, 白一色というイメージが強いので, カラーマップの作成の必要がなくなってしまうためでした. 最後まで読んでいただきありがとうございました.