residual.m

MDPSAS工具箱是马里兰大学开发的

function A = residual(A)

unpack(A)

b = 5.677e-8;   % Boltzmann constant, W/(m^2 K^4)
Ew = 0.7;      % wafer emissivity
Eq = 0.37;     % guard ring (quartz) emissivity

kcoeffw = [269 -0.585 3.75e-4];
kcoeffg = [1.53 -1.87e-3 5.8e-6];

Cpcoeffw = [307 1.54 -1.06e-3];
Cpcoeffg = [330 1.77 -9.45e-4];

kw = kcoeffw(1) + kcoeffw(2)*T + kcoeffw(3)*T^2;
Cpw = Cpcoeffw(1) + Cpcoeffw(2)*T + Cpcoeffw(3)*T^2;
dCpw = Cpcoeffw(2) + 2*Cpcoeffw(3)*T;

kg = kcoeffg(1) + kcoeffg(2)*T + kcoeffg(3)*T^2;
Cpg = Cpcoeffg(1) + Cpcoeffg(2)*T + Cpcoeffg(3)*T^2;
dCpg = Cpcoeffg(2) + 2*Cpcoeffg(3)*T;

k   = Wmask*(0.1*kw + 0.9*kg) + (1-Wmask)*kg;
k = bpsf(k,1);
Cp  = Wmask*(0.1*Cpw + 0.9*Cpg) + (1-Wmask)*Cpg;
dCp = Wmask*(0.1*dCpw + 0.9*dCpg) + (1-Wmask)*dCpg;

Cptot = Cp + dCp;
emiss = Wmask*Ew + (1-Wmask)*Eq;

Rp = delz*( (Dr*k)*(Dr*T) + k*(DDr*T) ) - b*emiss*T^4 + Ql + Qg + Qsh;

% B.C. residuals

Rp = setbval(Rp,Dr*T,'r','min');
Rp = setbval(Rp,k*(Dr*T)+b*emiss*T^4,'r','max');

A = set(A,'resid',Rp);