gusucode.com > 热力学传导matlab源码程序 > Code\CantileverHeatEqIterative.m

    
function T = CantileverHeatEqIterative(Ts, T0, Gts, AAu, ASi, APt, wSi, h, x);

% Thermal Properties
kSi = 16;
kAu = 317;
kPt = 71.6;


% Set up initial temperature matrix at ambient

for l = 1:length(x)

    T(l) = T0;
    
end

J = length(x);

dx = mean(diff(x));

T(1) = T0; % Sink Boundary Condition
TLast = ones(size(x));

for k = 1:J

    while abs( TLast(k) - T(k)) > 10e-10 

        for j = 2:J-1

            T(j) = [ 1/[ wSi(j)*h(j) + [2/(dx^2)]*[ kAu*AAu(j) + kSi*ASi(j) + kPt*APt(j) ]] ]*[ [ (T(j+1)+ T(j-1))/(dx^2) ]*[ kAu*AAu(j) + kSi*ASi(j) + kPt*APt(j)] + [ (T(j+1)-T(j-1))/(2*dx) ]*[ kAu*((AAu(j+1)-AAu(j-1))/(2*dx)) + kSi*((ASi(j+1)-ASi(j-1))/(2*dx)) + kPt*((APt(j+1)-APt(j-1))/(2*dx)) ] + wSi(j)*h(j)*Ts(j) ] ;

        end

        T(J) = [ Ts(J)*(Gts/dx) + (T(J-1)/(dx^2))*((kAu*AAu(J))+(kSi*ASi(J))+(kPt*APt(J))) ] / [ (Gts/dx) + (1/(dx^2))*((kAu*AAu(J))+(kSi*ASi(J))+(kPt*APt(J))) ]; % Neumann Boundary Condition at x = L controlled by Gts

        TLast = T;

        for i = 2:J-1

            T(i) = [ 1/[ wSi(i)*h(i) + [2/(dx^2)]*[ kAu*AAu(i) + kSi*ASi(i) + kPt*APt(i) ]] ]*[ [ (TLast(i+1)+ TLast(i-1))/(dx^2) ]*[ kAu*AAu(i) + kSi*ASi(i) + kPt*APt(i)] + [ (TLast(i+1)-T(i-1))/(2*dx) ]*[ kAu*((AAu(i+1)-AAu(i-1))/(2*dx)) + kSi*((ASi(i+1)-ASi(i-1))/(2*dx)) + kPt*((APt(i+1)-APt(i-1))/(2*dx)) ] + wSi(i)*h(i)*Ts(i) ] ;

        end

        T(J) = [ Ts(J)*(Gts/dx) + (TLast(J-1)/(dx^2))*((kAu*AAu(J))+(kSi*ASi(J))+(kPt*APt(J))) ] / [ (Gts/dx) + (1/(dx^2))*((kAu*AAu(J))+(kSi*ASi(J))+(kPt*APt(J))) ]; % Neumann Boundary Condition at x = L controlled by Gts

    end

end