熱モデルのGeome​tryの設定にimp​ortGeometr​y関数とmultis​phere関数とを使​うのでシミュレーショ​ンの挙動が異なるのは​なぜですか?

私はPartial Differential Equation Toolboxを使って二層の球からなるジオメトリを持つ熱モデルのシミュレーションをしています. 内側の球を内部熱源として機能させています.その中で, importGeometry関数でGeometryを設定した熱モデルとmultisphere関数でGeometryを設定した熱モデルとでシミュレーションの挙動が異なることを確認しました. 内側の球の境界条件が変化しているように思えます. multisphere関数で設定した方は, 内側の球の熱をそのまま伝えているのに対し, importGeometry関数で設定した方は, 内側の球の境界で断熱が起きているような挙動になります. どうしてこのようなことが起きるのでしょうか.

%Using multisphere function
%Create thermalmodel
thermalmodel = createpde('thermal','transient');
gm = multisphere([3.15*10^-3,1*10^-2]);
thermalmodel.Geometry = gm;
%Visualize geometry
figure
subplot(1,2,1)
pdegplot(thermalmodel,'FaceAlpha',0.25,'CellLabel','on')
title('Geometry with Cell Labels')
subplot(1,2,2)
pdegplot(thermalmodel,'FaceAlpha',0.25,'FaceLabel','on')
title('Geometry with Face Labels')
%Set physical parameter
generateMesh(thermalmodel,'Hmax',1*10^-3);
thermalProperties(thermalmodel,'Cell',1,'ThermalConductivity',0.778,'MassDensity',1.66*10^3,'SpecificHeat',2.54*10^3);
thermalProperties(thermalmodel,'Cell',2,'ThermalConductivity',0.642,'MassDensity',1*10^3,'SpecificHeat',3.72*10^3);
internalHeatSource(thermalmodel,6.15*10^6,'cell',1);
thermalBC(thermalmodel,'Face',2,'ConvectionCoefficient',25,'AmbientTemperature',25);
thermalIC(thermalmodel,25);
%Calculate PDE
tlist=0:1:300;
result = solve(thermalmodel,tlist);
Tmin=0;
Tmax=max(result.Temperature(:,end));
%Visualize by contour
[YG,ZG] = meshgrid(linspace(-1*10^-2,1*10^-2,100),linspace(-1*10^-2,1*10^-2,100));
XG = zeros(size(YG));
tIndex = [6,51,101,151,201,301];
varNames = {'Time_index','Time_step'};
index_step = table(tIndex.',tlist(tIndex).','VariableNames',varNames);
disp(index_step);
TG = interpolateTemperature(result,XG,YG,ZG,tIndex);
t = linspace(0,2*pi);
ylayer1 = 3.15*10^-3*cos(t); zlayer1 = 3.15*10^-3*sin(t);
ylayer2 = 1*10^-2*cos(t); zlayer2 = 1*10^-2*sin(t);
figure('Position',[10,10,1000,550]);
for i = 1:numel(tIndex)
    subplot(2,3,i)
    contour(YG,ZG,reshape(TG(:,i),size(YG)),'ShowText','on')
    colormap turbo 
    c=colorbar;
    c.Label.String = 'Temperature(℃)';
    title(['Temperature at Time ' num2str(tlist(tIndex(i)))]);
    hold on
    caxis([Tmin,Tmax])
    axis equal
    plot(ylayer1,zlayer1,'k','LineWidth',1.5)
    plot(ylayer2,zlayer2,'k','LineWidth',1.5)
   end

↓Result

%Using importGeometry function
%Create thermalmodel
thermalmodel = createpde('thermal','transient');
gm = importGeometry('ball_simu.stl');
thermalmodel.Geometry = gm;
%Visualize geometry
figure
subplot(1,2,1)
pdegplot(thermalmodel,'FaceAlpha',0.25,'CellLabel','on')
title('Geometry with Cell Labels')
subplot(1,2,2)
pdegplot(thermalmodel,'FaceAlpha',0.25,'FaceLabel','on')
title('Geometry with Face Labels')
%Set physical parameter
generateMesh(thermalmodel,'Hmax',1*10^-3);
thermalProperties(thermalmodel,'Cell',2,'ThermalConductivity',0.778,'MassDensity',1.66*10^3,'SpecificHeat',2.54*10^3);
thermalProperties(thermalmodel,'Cell',1,'ThermalConductivity',0.642,'MassDensity',1*10^3,'SpecificHeat',3.72*10^3);
internalHeatSource(thermalmodel,6.15*10^6,'cell',2);
thermalBC(thermalmodel,'Face',1,'ConvectionCoefficient',25,'AmbientTemperature',25);
thermalIC(thermalmodel,25);
%Calculate PDE
tlist=0:1:300;
result = solve(thermalmodel,tlist);
Tmin=0;
Tmax=max(result.Temperature(:,end));
%Visualize by contour
[YG,ZG] = meshgrid(linspace(-1*10^-2,1*10^-2,100),linspace(-1*10^-2,1*10^-2,100));
XG = zeros(size(YG));
tIndex = [6,51,101,151,201,301];
varNames = {'Time_index','Time_step'};
index_step = table(tIndex.',tlist(tIndex).','VariableNames',varNames);
disp(index_step);
TG = interpolateTemperature(result,XG,YG,ZG,tIndex);
t = linspace(0,2*pi);
ylayer1 = 3.15*10^-3*cos(t); zlayer1 = 3.15*10^-3*sin(t);
ylayer2 = 1*10^-2*cos(t); zlayer2 = 1*10^-2*sin(t);
figure('Position',[10,10,1000,550]);
for i = 1:numel(tIndex)
    subplot(2,3,i)
    contour(YG,ZG,reshape(TG(:,i),size(YG)),'ShowText','on')
    colormap turbo 
    c=colorbar;
    c.Label.String = 'Temperature(℃)';
    title(['Temperature at Time ' num2str(tlist(tIndex(i)))]);
    hold on
    caxis([Tmin,Tmax])
    axis equal
    plot(ylayer1,zlayer1,'k','LineWidth',1.5)
    plot(ylayer2,zlayer2,'k','LineWidth',1.5)
   end

↓Result

Answer 1

[Hiroyuki Hishida]

佐藤様、

generateMeshで生成される形状が異なるためです。例えばthermalproperty設定の前あたりに、以下を差し込んで実行してみてください。すると、multiSphereとimportのときで、表示される形状（点群）が異なることが確認できると思います。

myMesh = generateMesh(thermalmodel,'Hmax',0.01);
figure
pdemesh(myMesh,'FaceAlpha',0)
title('generateMeshで生成される形状')

佐藤様のスクリプトでもそうですが、STLをimportしてもそれをそのまま解析に渡すことができず、generatemeshします。このgenerateMeshは、与えられた形状の主に頂点情報から有限要素を生成します。そのため、与える形状情報はstlでなくても、csg表現でも、alpha shapeにも対応します。

この与える形状情報ですが、multisphereの場合、球の中に球があるカタチにはなっておらず、あくまでモノとしての形状は一つの球です。一方、佐藤様のスクリプトにおいては元のSTLにおいても球の中に球を作成されています。ここが違いとしてあります。

シンプルな例題としては以下が参考になるかと思います。

https://jp.mathworks.com/help/pde/ug/deflection-analysis-of-a-bracket.html

では単一の球をSTLで与えて内部に球をCELLで追加したい場合どうするかですが、例えば以下はいかがでしょうか。

https://jp.mathworks.com/help/pde/ug/discretegeometry.addcell.html

どうぞよろしくお願いします。

菱田

[和奈 佐藤]

参考にさせて頂いたところ, うまくいきました.

ありがとうございます!

[Hiroyuki Hishida]

うまくいって良かったです。

菱田