スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

FEM:2次元定常熱伝導解析

久々にお勉強です。
Personal CAE Project
こちらのサイトの熱伝導解析の例題1をMatlabでやってみました。


FEMthermal20140130.png
まずはサイトに大感謝です。
熱伝導解析に関しての知識はほぼゼロでしたが
ごく短時間で例題をプログラミングするところまで辿り着けました。
その精神を見習って私もスクリプトを公開しておきます。
(雑なコードですが結果は出ますのでご勘弁を)

ここから下です。

clear all
close all

%--- Input Value --- %
Lx = 100;
Ly = 50;
d_x = 50;
d_y = 25;

ramda = 0.225;

%--- Boundary Contidion ---%
BC_node = [];
BC_T = [];

bc1 = 1:d_x;
bc1_T = zeros(1,length(bc1));
bc2 = d_x+2:d_x+1:(d_x+1)*d_y;
bc2_T = zeros(1,length(bc2)) + 50;
bc3 = d_x+1:d_x+1:(d_x+1)*d_y;
bc3_T = zeros(1,length(bc3)) + 100;
bc4 = (d_x+1)*d_y+1:(d_x+1)*(d_y+1);
bc4_T = zeros(1,length(bc4)) + 100;

BC_node = [BC_node bc1 bc2 bc3 bc4 ];
BC_T = [BC_T bc1_T bc2_T bc3_T bc4_T ];

%--- Node coordinate ---%
N_node = (d_x+1) * (d_y+1);

Coord_node = zeros(N_node, 2);
Coord_node(1:d_x+1,1) = [0:d_x]*(Lx/d_x);
for i = 1:d_y
Coord_node([1:d_x+1]+(d_x+1)*i,1) = Coord_node(1:d_x+1,1);
Coord_node([1:d_x+1]+(d_x+1)*i,2) = Ly/d_y * i;
end

%--- Element ---%
N_ele = d_x * d_y * 2;
Node_ele = zeros(N_ele, 3);
Node_ele(1,:) = [1 2 (d_x+1)+1];
Node_ele(2,:) = [2 (d_x+1)+2 (d_x+1)+1];
for i = 2:d_x
Node_ele(1+(i-1)*2,:) = Node_ele(1,:) + (i-1);
Node_ele(2+(i-1)*2,:) = Node_ele(2,:) + (i-1);
end
for i = 2:d_y
Node_ele([1:d_x*2]+(i-1)*d_x*2,:) = Node_ele([1:d_x*2],:) + (d_x+1)*(i-1);
end


%--- [K] matrix ---%
K = zeros(N_node);
for e = 1:N_ele
Ke = zeros(3);
n1 = Node_ele(e,1);
n2 = Node_ele(e,2);
n3 = Node_ele(e,3);

x1 = Coord_node(n1,1);
x2 = Coord_node(n2,1);
x3 = Coord_node(n3,1);
y1 = Coord_node(n1,2);
y2 = Coord_node(n2,2);
y3 = Coord_node(n3,2);

b1 = y2 - y3;
b2 = y3 - y1;
b3 = y1 - y2;
c1 = x2 - x3;
c2 = x3 - x1;
c3 = x1 - x2;

delta_e = 0.5*( (x2-x1)*(y3-y1)...
-(x3-x1)*(y2-y1) );

Ke(1,1) = b1*b1 + c1*c1;
Ke(1,2) = b1*b2 + c1*c2;
Ke(1,3) = b1*b3 + c1*c3;
Ke(2,2) = b2*b2 + c2*c2;
Ke(2,3) = b2*b3 + c2*c3;
Ke(3,3) = b3*b3 + c3*c3;
Ke(2,1) = Ke(1,2);
Ke(3,1) = Ke(1,3);
Ke(3,2) = Ke(2,3);

Ke = ramda * 0.25 / delta_e * Ke;

K([n1 n2 n3], [n1 n2 n3]) = K([n1 n2 n3], [n1 n2 n3]) + Ke;

end

%--- Solove [K]*{T} = {0} ---%
Free_node =[];
for i = 1:N_node
if min( abs(BC_node - i) ) == 0
else
Free_node = [Free_node i];
end
end

b = - K(:,BC_node)*BC_T.';
K2 = K(:,Free_node);

T = zeros(N_node,1);
T(BC_node,1) = BC_T.';
T(Free_node,1) = K2\b;

figure(1)
for e = 1:N_ele
n1 = Node_ele(e,1);
n2 = Node_ele(e,2);
n3 = Node_ele(e,3);

x(1) = Coord_node(n1,1);
x(2) = Coord_node(n2,1);
x(3) = Coord_node(n3,1);
y(1) = Coord_node(n1,2);
y(2) = Coord_node(n2,2);
y(3) = Coord_node(n3,2);

patch(x,y,T([n1 n2 n3]) );
end
colorbar
set(gca,'DataAspectRatio',[1 1 1])
スポンサーサイト

コメントの投稿

非公開コメント

カレンダー
09 | 2017/10 | 11
1 2 3 4 5 6 7
8 9 10 11 12 13 14
15 16 17 18 19 20 21
22 23 24 25 26 27 28
29 30 31 - - - -
プロフィール

uglab

Author:uglab
ソフト開発会社

【Twitter】
フォローお願いします

【You tube動画一覧】
チャンネル登録お願いします


最新記事
最新コメント
最新トラックバック
月別アーカイブ
カテゴリ
FC2カウンター
検索フォーム
スポンサードリンク
    【スポンサードリンク】

リンク
RSSリンクの表示
ブロとも申請フォーム

この人とブロともになる

QRコード
QRコード
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。