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