cc1 = r*dr/dt; viscD = 2*eta.*(((uf(2:end) - uf(1:end-1))/dr).^2 + (u./r).^2 - 1/3*((uf(2:end).*rf(2:end) - uf(1:end-1).*rf(1:end-1))/dr./r).^2); AT = spdiags([-lambdaf(2:nF).*rf(2:nF)/dr - max(mFlux(1:end-1),0).*[0;cp(1:end-1)], (cp.*cc1.*rho + [lambdaf(2:nF-1).*rf(2:nF-1)/dr;0] + lambdaf(1:nF-1).*rf(1:nF-1)/dr) + max(mFlux(2:end),0).*cp + max(-mFlux(1:end-1),0).*cp , -lambdaf(1:nF-1).*rf(1:nF-1)/dr- max(-mFlux(1:end-1).*cp,0)],-1:1,nC,nC); %AT = spdiags([-lambdaf(2:end).*rf(2:end)/dr - max(mFlux(1:end-1).*[0;cp(1:end-1)]), (cp.*cc1.*rho + [lambdaf(2:end-1).*rf(2:end-1)/dr;lambdaf(end)*rf(end)/0.5/dr] + lambdaf(1:end-1).*rf(1:end-1)/dr) + max(mFlux(2:end),0).*cp + max(-mFlux(1:end-1),0).*cp, -lambdaf(1:end-1).*rf(1:end-1)/dr - max(-mFlux(1:end-1).*cp,0)],-1:1,nC,nC); %T = 400 end %beta = betaFunc(T,p); beta = -1; bT = cc1.*cp0.*T0.*rho0 + rho.*T.*r*dr.*((cp - cp0)/dt + u.*(cpf(2:end) - cpf(1:end-1))/dr) - beta.*((p - p0)/dt+ u.*gradp).*r*dr + viscD.*r*dr; %bT(end) -= max(-mFlux(end),0)*hu; %bT(end) += lambdaf(end)*rf(end)/0.5/dr*400; bT += initSource.*r*dr;% + max(0,polyval(EgPol,t))*tanhProfile'.*r*dr + Qc.*r*dr/dt.*rho.*(1 - burnt).*(T>Tadb); Ty = AT\bT; TInitRes = res_compute(AT,bT,T); T = (1 - URFt)*T + URFt*Ty; TFinalRes = res_compute(AT,bT,T); if cc == 130 URFp = 0.8; URFt = 0.5; elseif cc == 80 URFp = 0.1; URFt = 0.1; %URFi_u = 0.2; elseif cc == 50 URFp = 0.3; URFt = 0.3; %URFi_u = 0.4; elseif cc >=170 uInitRes TInitRes pInitRes keyboard endif [cp,lambda,M,rho,eta] = updtProp(T,p); %T = T*(1 - URFt) + URFt*TFromh(T,h,p); %T = TFromh(T,h,p); rhof = interpolate(rho); cpf = interpolate(cp); %eta = 0; %T = h./cp; %cp = 1000;% + T*.001; %lambda = 0.02 + T.*0;% + T*.0001; %M = 29 + T.*0;