%clear all recov = 0; debug = 0; dbf = 15; itmax = 200; tol = 1e-3; icri = 0; dg = 2e-3; pInit = 101325; Tu = 298; gamma = 1.37; Cbd = 200; Ebd = (4.3 + pInit/1e5/Tu*(136 + 324*dg*1000))^2/Cbd^2/dg/1000; Ti = Tu*(1 + (60000/Tu - 1)/gamma); di = 2*sqrt((gamma - 1)/gamma*Ebd/pInit/dg/(1 - Tu/Ti)/pi); ri = 0.5*di; Tadb = 1892; Qc = 1934670; Le = 1.7; Su0 = 0.203; dr = 30e-6; rmin = 0;%dr; rmax = 0.24e-2; dt = 2e-11;%0.2e-9; nC = round((rmax - rmin)/dr); nF = nC + 1; t_end = 100e-6; %rk = zeros(ceil(t_end/dt)+1,1); rk(1) = ri; rp = rk; rp2 = rp; rpCor = rp2; Energy = 54.8; switch Energy case 2 EgPol = [-1.9231e+11 3.4208e+07 -4.0267e+04 1.2831e+01]; case 11.5 EgPol = [-1.7604e+10 2.6177e+07 -4.0369e+04 2.7137e+01]; case 54.8 EgPol = [-2.9110e+08 1.7790e+06 -9.9874e+03 2.9136e+01]; case 12.4 EgPol = [-9.6989e+09 1.8887e+07 -3.9422e+04 2.8194e+01]; endswitch %Egdot = max(0,polyval(EgPol,t)); %%electEnergy eta_e = 0.5; rf = (rmin:dr:rmax)'; r = 0.5*(rf(1:end-1)+rf(2:end)); %rf = rf(2:end); %Solve diffusion equation from post-breakdown to T = 3 x Tadb uf0 = ones(nF,1)*1e-15*0; u = uf0(2:end); uf = uf0; ufavg0 = uf0; ufavg = ufavg0; u0 = u; cp = 1000*ones(nC,1); cp0 = cp; pf = ones(nF,1).*pInit; p = pf(2:end); pp = p*0; p0 = p; pf = [p(1);0.5*(p(1:end-1)+p(2:end));pInit];%p(end-1) + (p(end) - p(end-1))*1.5]; gradp = (pf(2:end) - pf(1:end-1))/dr; gradpf = [0;(p(2:end)-p(1:end-1))/dr;(pInit-p(end))/.5/dr]; gradpfavg = [gradp(1);0.5*(gradp(1:end-1)+gradp(2:end));gradp(end)]; %uf0 = 0.9*uf0 + 0.1*uf0y; %plot(uf0) % %(qf - [0;qf(1:end-1)] + uFlux.*uf0 - [0;uFlux(1:end-1).*uf0(1:end-1)]); %keyboard eta = 1e-5*ones(nC,1); M = ones(nC,1)*29; lambda = 0.025*ones(nC,1); lambdau = 0.025; T = ones(80,1)*300; rho = p/8314*29./T; rho0 = rho; rhof = interpolate(rho); etaf = interpolate(eta); Mf = interpolate(M); rhof0 = rhof; lambdaf = [lambda(1);0.5*(lambda(1:end-1)+lambda(2:end));lambdau]; M0 = M; uf = uf0; mFlux = rf.*uf.*rhof; cpf = [cp(1);0.5*(cp(1:end-1)+cp(2:end));cp(end)]; t = 0; i = 0; subplot (3, 1, 1); hold off plot(0) subplot (3, 1, 2) hold off plot(0) subplot (3, 1, 3) hold off plot(0) h = 1000*ones(nC,1); uDel = 10 + u*0; pDel = uDel; hDel = pDel; if recov load -binary "recover.m" t h u uf ufavg T rho p M i h0 = h; u0 = u; uf0 = uf; T0 = T; rho0 = rho; p0 = p; M0 = M; endif pInitRes = 0; uInitRes = 0; %rho = p/8314*Mu/Tu; rho0 = rho; p = rho*8314./M.*T; p0 = p; %h(1:10) *= 10; h0 = h; uInitRes = 0; pInitRes = 0; dt = 1e-7; while t= 100000 uInitRes TInitRes pInitRes uy(1) uf(1) hy(1) mFlux(1) ppy(1) keyboard endif u_1 = u; p_1 = p; h_1 = h; %uEqn %eta = 0; momentumEqStaggered %energyEquation %TInitRes %keyboard % continuityEquation % pressureCorrEq %enthalpyEquation %keyboard TEquation %pEqOF %pCorrEqn %p = p0.*M0./T0.*T./M; %keyboard if 0 Au = spdiags(ones(nC,1),0,nC,nC) - spdiags(ones(nC-1,1),-1,nC,nC); bu = -(rho - rho0)/dt*dr.*r; mFlux(2:nF) = Au\bu; %mFlux = min(mFlux,.02); u = 0.5*(mFlux(1:end-1)+mFlux(2:end))./rho./r; uf = mFlux./rhof./[1;rf(2:end)]; endif if 0 uDel_1 = uDel; uDel = u - u_1; pDel_1 = pDel; pDel = p - p_1; hDel_1 = hDel; hDel = h - h_1; uLamb = norm(uDel)/norm(uDel_1 +1e-10); pLamb = norm(pDel)/norm(pDel_1 +1e-10); hLamb = norm(hDel)/norm(hDel_1 +1e-10); uErr = norm(uDel)/(uLamb-1); pErr = norm(pDel)/(pLamb-1); hErr = norm(hDel)/(hLamb-1); Err = max(abs([uErr,pErr,hErr])); endif Err = max([uInitRes;pInitRes;TInitRes]); %if uInitRes < 1e-4 && max(abs(pp./p)) < 1e-2 && TInitRes < 1e-5 if Err < tol % || max(abs(pp))<1 out = 1; endif cc++; %keyboard until ((cc >= itmax) || (out == 1)) % || max(abs(pp))<1e-3 %keyboard %plot (r,T) % pause(0.5) % keyboard %uInitRes %TInitRes %pInitRes % (rho-rho0)/dt.*r*dr + mFlux(2:end) - mFlux(1:end-1); % (rho.*h - rho0.*h0) if i == 100 dt = 5e-10; elseif i == 150 dt = 2e-9; elseif i == 200 dt = 6e-9; endif if cc == itmax disp("Falla de convergencia") keyboard endif if 0 if i > 100 cri = max([max(abs(T-T0)./T0),max(abs(p-p0)./p0),max(abs(u-u0)./(u0 + 10))]); %cri = max(abs(T-T0)./T0); if cri > .3 && i > icri + 5 icri = i; %keyboard dt /= 1.5; elseif cri < .05 && i > icri + 5 icri = i; %keyboard dt *= 1.1; endif endif dt = min(max(dt,1e-12),20e-6); endif if 0 if i == 50 || i ==100 || i == 150 || i == 300 || i == 450 || i == 600 || i == 750 || i == 900 dt *= 5; endif if min(p) < pInit*.7 dt = 1e-9; endif endif %if p(1) < 295000 && t>100*dt % keyboard % endif h0 = h; u0 = u; uf0 = uf; ufavg0 = ufavg; T0 = T; rho0 = rho; p0 = p; M0 = M; cp0 = cp; save -binary "recover.m" t h u uf ufavg T rho p M i if (mod(i,1) == 0) && 1 subplot (3, 1, 1); plot (r,T); title('T') subplot (3, 1, 2) plot (r,u); title('u') subplot (3, 1, 3) plot (r,p); title('p') pause(0.5) endif % keyboard if (mod(i,1) == 0)&&0 subplot (3, 1, 1); plot (r(1:20),T(1:20)); title('T') subplot (3, 1, 2) hold on plot (rf(1:20),uf(1:20)); title('u') subplot (3, 1, 3) plot (r(1:20),p(1:20)); title('p') pause(0.5) endif burnt = max(burnt,T>=Tadb); if 0 Tfil = sum(T>Tadb); rk(i+1) = r(Tfil) + (Tadb - T(Tfil))/(T(Tfil+1)-T(Tfil))*(r(Tfil+1) - r(Tfil)); sk = (rk(i+1) - rk(i))/dt; if (max(T) > 3000) eta_e = 0.36 + 0.14*sk^3/(sk^3 + 500); else eta_e = 0.08 + 0.22*sk^3/(sk^3 + 700); endif %D = %Y = Le*rho*cp/lambda RR = 2*1.68e-3/rpCor(i); RRR = 1 + 1/RR + 1/RR^2 +2/3/RR^3; rp(i+1) = rp(i) + dt*(sk + 6.3*Su0*max(0.7,rp(i)/(rp(i)+2*1.68e-3))); rp2(i+1) = rp2(i) + dt*(sk + 6.3*Su0*max(0.9,rp(i)/(rp(i)+2*1.68e-3))); rpCor(i+1) = rpCor(i) + dt*(max(sk, 6.3*Su0/(RRR))); endif endwhile if 0 subplot (3, 2, 1); plot (T(1:10)); title('T') subplot (3, 2, 2) plot (u(1:10)); title('u') subplot (3, 2, 3) plot (mFlux(1:10)); title('mFlux') subplot (3, 2, 4) plot (p(1:10)); title('p') subplot (3, 2, 5) plot (pp(1:10)); title('pp') subplot (3, 2, 6) plot (gradp(1:10)); title('gradp') %endif subplot (1, 1, 1) plot (r,u); hold on plot (rf,uf); hold off endif