有限元方程：

其中

具体计算公式

对于二维三节点三角形的情况，各公式计算如下

ij边界上

jm边界上

mi边界上

具体计算同b)，但是由于这里是稳定温度场，故可以不必考虑。

例. 如图所示的正方形平板，确定其温度分布。已知λ=10 W/(m.℃)，对流换热系数，自由流动温度，对流发生在4-5表面，1-5、2-4表面绝热，1-2表面温度保持100℃，板厚10mm。

调用程序，运行结果如下

第1个单元:

Ke=

0.05 0 -0.05

0 0.05 -0.05

-0.05 -0.05 0.1

Pe=

0

0

0

第2个单元:

Ke=

0.05 0 -0.05

0 0.05 -0.05

-0.05 -0.05 0.1

Pe=

0

0

0

第3个单元:

Ke=

0.183333333333333 0.0666666666666667 -0.05

0.0666666666666667 0.183333333333333 -0.05

-0.05 -0.05 0.1

Pe=

10

10

0

第4个单元:

Ke=

0.05 0 -0.05

0 0.05 -0.05

-0.05 -0.05 0.1

Pe=

0

0

0

[K]*T={P}

其中

[K]=

0.1 0 -0.1 0 0

0 0.1 -0.1 0 0

-0.1 -0.1 0.4 -0.1 -0.1

0 0 -0.1 0.233333333333333 0.0666666666666667

0 0 -0.1 0.0666666666666667 0.233333333333333

{P}=

0

0

0

10

10

各节点温度为:

100

100

80

60

60

 

具体程序如下

main.m文件

clc
clear all
close all

switch 2
    case 1 % 书上例题
        verts = [4 6;
            12 8;
            8 10;];
        eleindex = [1 2 3;];
        edgeindex = [1 2;
            2 3;]; % 只标记有热量交换的边
        alpha = [15 0 40;
            10 0 40;]; % 表面对流换热系数alpha、表面热流q、温度tf
        lambda = 60;
        t = 1;
        qv = 50;
    case 2 % 作业第4题
        verts = [0 2;
            0 0;
            1 1;
            2 0;
            2 2;];
        eleindex = [1 2 3;
            2 4 3;
            4 5 3;
            5 1 3;];
        edgeindex = [1 2;
            4 5;];
        alpha = [0 0 100;
            20 0 50;]; % 表面对流换热系数alpha、表面热流q、温度tf
        lambda = 10;
        t = 0.01;
        qv = 0;
end

% 构造FEM
hFEM.vertcount = size(verts, 1);
hFEM.verts = verts;
hFEM.elecount = size(eleindex, 1);
hFEM.eleindex = eleindex;
hFEM.edgeindex = edgeindex;
hFEM.alpha = alpha;
hFEM.lambda = lambda;
hFEM.t = t;
hFEM.qv = qv;

hFEM = FEM(hFEM);
for i = 1:hFEM.elecount
    fprintf('第%d个单元:\n', i);
    disp('Ke=')
    disp(hFEM.eletruct(i).K)
    disp('Pe=')
    disp(hFEM.eletruct(i).P)
end
fprintf('\n\n\n');

% 求解
fprintf('[K]*T={P}\n其中\n');
K = hFEM.K;
P = hFEM.P;
disp('[K]=')
disp(K)
disp('{P}=')
disp(P)
fprintf('\n\n\n');

fixind = [1 2];
unfixind = [3 4 5];
K11 = K(fixind, fixind);
K12 = K(fixind, unfixind);
K21 = K(unfixind, fixind);
K22 = K(unfixind, unfixind);
T1 = [100; 100];
T2 = K22\(P(unfixind)-K21*T1);
T(fixind) = T1;
T(unfixind) = T2;
T = T';
fprintf('各节点温度为:\n');
disp(T)

 

FEM.m文件

% 计算各单元的刚度矩阵
function hFEM = FEM(hFEM)
% 处理每个单元
for i = 1:hFEM.elecount
    hFEM.eletruct(i) = FEMele(hFEM, i);
end

% 计算合成刚度矩阵
hFEM.K = FEMCalcCombiningK(hFEM);

% 计算热载荷
hFEM.P = zeros(hFEM.vertcount, 1);
for i = 1:hFEM.elecount
    [hFEM.eletruct(i).P hFEM.eletruct(i).P1 hFEM.eletruct(i).P2 hFEM.eletruct(i).P3] = CalcEleThermalLoad(hFEM, i);
    hFEM.P(hFEM.eleindex(i, :)) = hFEM.P(hFEM.eleindex(i, :)) + hFEM.eletruct(i).P;
end
end

function hFEMele = FEMele(hFEM, ind)
hFEMele.ind = hFEM.eleindex(ind, :);
hFEMele.ele = hFEM.verts(hFEM.eleindex(ind, :), :);
hFEMele.B = FEMCalcB(hFEMele);
hFEMele.D = FEMCalcD(hFEM);
hFEMele.K1 = FEMCalcK1(hFEM, hFEMele);
hFEMele.K2 = FEMCalcK2(hFEM, hFEMele);
hFEMele.K = hFEMele.K1 + hFEMele.K2;
end

% 计算几何矩阵
function B = FEMCalcB(hFEMele)
ele = hFEMele.ele;
A(:, 2:3) = ele;
A(:, 1) = 1;
invA = inv(A);

B = invA(2:3, :);
end

% 计算应力矩阵
function D = FEMCalcD(hFEM)
D = hFEM.lambda*eye(2);
end

% 计算单元刚度矩阵K1
function K1 = FEMCalcK1(hFEM, hFEMele)
triarea = CalcTriangleArea(hFEMele.ele);

K1 = hFEM.t*triarea*hFEMele.B'*hFEMele.D*hFEMele.B;
end

% 计算单元刚度矩阵K2
function K2 = FEMCalcK2(hFEM, hFEMele)
lij = norm(hFEMele.ele(2, :) - hFEMele.ele(1, :));
Sij = hFEM.t*lij;
aij = FindAlpha(hFEMele.ind([1 2]), hFEM);
K2ij = Sij*aij/6*[2 1 0; ...
    1 2 0; ...
    0 0 0;];

ljm = norm(hFEMele.ele(3, :) - hFEMele.ele(2, :));
Sjm = hFEM.t*ljm;
ajm = FindAlpha(hFEMele.ind([2 3]), hFEM);
K2jm = Sjm*ajm/6*[0 0 0; ...
    0 2 1; ...
    0 1 2;];

lmi = norm(hFEMele.ele(3, :) - hFEMele.ele(1, :));
Smi = hFEM.t*lmi;
ami = FindAlpha(hFEMele.ind([3 1]), hFEM);
K2mi = Smi*ami/6*[2 0 1; ...
    0 0 0; ...
    1 0 2;];

K2 = K2ij + K2jm + K2mi;
end

% 找指定边的对流换热系数
function a = FindAlpha(ind, hFEM)
for i = 1:length(hFEM.edgeindex)
    if all(ind == hFEM.edgeindex(i, :)) || all(ind([2, 1]) == hFEM.edgeindex(i, :))
        a = hFEM.alpha(i, 1);
        return;
    end
end
a = 0;
end

% 找指定边的表面热流
function q = FindQ(ind, hFEM)
for i = 1:length(hFEM.edgeindex)
    if all(ind == hFEM.edgeindex(i, :)) || all(ind([2, 1]) == hFEM.edgeindex(i, :))
        q = hFEM.alpha(i, 2);
        return;
    end
end
q = 0;
end

% 找指定边的温度
function tf = FindTf(ind, hFEM)
for i = 1:length(hFEM.edgeindex)
    if all(ind == hFEM.edgeindex(i, :)) || all(ind([2, 1]) == hFEM.edgeindex(i, :))
        tf = hFEM.alpha(i, 3);
        return;
    end
end
tf = 0;
end

% 计算三角形面积
function area = CalcTriangleArea(ele)
A(:, 2:3) = ele;
A(:, 1) = 1;
area = 0.5*det(A);
end

% 计算合成矩阵
function K = FEMCalcCombiningK(hFEM)
vertcount = hFEM.vertcount;
elecount = hFEM.elecount;
eleindex = hFEM.eleindex;

eleK = zeros(3, 3, elecount);
for k = 1:elecount
    eleK(:, :, k) = hFEM.eletruct(k).K;
end

K = zeros(vertcount);
for i = 1:vertcount
    for j = 1:vertcount
        for k = 1:elecount
            indi = find(eleindex(k, :)==i);
            indj = find(eleindex(k, :)==j);
            if ~isempty(indi) && ~isempty(indj)
                K(i, j) = K(i, j) + eleK(indi:indi, indj:indj, k);
            end
        end
    end
end
end

% 计算各单元的热载荷
function [P P1 P2 P3] = CalcEleThermalLoad(hFEM, ind)
hFEMele = hFEM.eletruct(ind);
aij = FindAlpha(hFEMele.ind([1 2]), hFEM);
qij = FindQ(hFEMele.ind([1 2]), hFEM);
tfij = FindTf(hFEMele.ind([1 2]), hFEM);
ajm = FindAlpha(hFEMele.ind([2 3]), hFEM);
qjm = FindQ(hFEMele.ind([2 3]), hFEM);
tfjm = FindTf(hFEMele.ind([2 3]), hFEM);
ami = FindAlpha(hFEMele.ind([3 1]), hFEM);
qmi = FindQ(hFEMele.ind([3 1]), hFEM);
tfmi = FindTf(hFEMele.ind([3 1]), hFEM);
lij = norm(hFEMele.ele(2, :) - hFEMele.ele(1, :));
Sij = hFEM.t*lij;
ljm = norm(hFEMele.ele(3, :) - hFEMele.ele(2, :));
Sjm = hFEM.t*ljm;
lmi = norm(hFEMele.ele(3, :) - hFEMele.ele(1, :));
Smi = hFEM.t*lmi;

P1 = hFEM.qv*CalcTriangleArea(hFEMele.ele)/3*[1; 1; 1];
P2 = qij*Sij/2*[1; 1; 0] + qjm*Sjm/2*[0; 1; 1] + qmi*Smi/2*[1; 0; 1];
P3 = tfij*aij*Sij/2*[1; 1; 0] + tfjm*ajm*Sjm/2*[0; 1; 1] + tfmi*ami*Smi/2*[1; 0; 1];

P = P1 - P2 + P3;
end