Решение ошибок с использованием pdepe для решения уравнения теплопроводности в осесимметричных полярных координатах в многослойном кольце

Question

Я пытаюсь смоделировать теплопроводность в витке медного провода. Я считаю, что функция pdepe является подходящей, поскольку проблема заключается в уравнении вынужденного 1D тепла в цилиндрических полярных координатах.

Как рекомендовано в другом месте на этом форуме, я моделирую весь составной домен (изоляция-медь-изоляция) как один домен с переменными свойствами в пространстве.

Уравнение в синтаксисе pdepe: (1 / alpha (r))du / dt = (1 / r)(r * du / dr) + q (r) / k(r)

Температура u(r,t), коэффициент диффузии альфа (r), объемная скорость тепловыделения q(r), проводимость k(r)

Он имеет внешние граничные условия du/dr = 0 (изолированный) на внутренней поверхности и k(r)du / dr + h(u-u0) = 0 (конвективный) на внешней поверхности.

Начальные условия везде нулевые.

Я сравнил вывод как t -> inf с аналитическим стационарным решением (прилагается), которое совершенно другое. Повышение температуры должно быть порядка 120 градусов С в устойчивом состоянии.

Устойчивое состояние профиля решения

Выходы из кода (соответствующие выдержки прилагаются ниже), однако, все O (1). Я был бы признателен, если бы кто-нибудь мог пролить свет на то, что я делаю неправильно, или если pdepe просто не подходит здесь. Заранее большое спасибо.

 function [sol,systInfo]=ThreeLayerHeat
   %uses pdepe()

  %% geometry

  rInner=0.05;
  rOuter=0.1;
  tIns=50e-6; %insulation thickness

  r2=rInner+tIns; %internal boundaries
  r3=rOuter-tIns;
  tDomain=rOuter-rInner;
  tCond=tDomain-tIns; %one conductor layer
  m = 1; %CPCs

  %% material properties

  %copper
  k_cond=394;
  rho_cond=8960;
  cp_cond=390;

  %mica
  k_ins=0.5;
  rho_ins=2800;
  cp_ins=880;

  alpha_cond=k_cond/(rho_cond*cp_cond);
  alpha_ins=k_ins/(rho_ins*cp_ins);

  %%Internal generation

  Jlimit=4e6; %4A/mm^2
  ACfactor=1.2; %CF design of electrical machines
  rho_conductor=1.75e-8; %electrical resistivity Table 9.3 293K %at reference level

  %%boundary details

  Tinf=0; %fluid temp on outer boundary
  h2=100; %HTC on outer boundary

  %%Characteristic (Fo) times

  tau_chic_cond=tCond^2/alpha_cond;
  tau_chic_ins=tIns^2/alpha_ins;
  tau_chic_limiting=max(tau_chic_cond,tau_chic_ins);

  %%discretisation

  %spatial
  npoints=50;

  x=linspace(rInner,rOuter,npoints);
  %ensure mesh points exist at internal boundaries r2 and r3
  x=unique(sort(horzcat(x,r2,r3)));

  %temporal
  nt=10;
  nPeriods=50; %how many times limiting (longest) tau (>1 to reach equilibrium)

  %linspace
  t0=0;
  tf=t0+tau_chic_limiting*nPeriods;
  t=linspace(t0,tf,nt);

  if t0~=0
      t=horzcat(0,t); %add 0 padding if doesn't start at 0 to allow initialisation
  end

  %logspace - uncomment this to use log time
  % logi=-2; %initial point
  % logf=log10(tf);
  % t=logspace(logi,logf,nt);
  %  t=horzcat(0,t); 

    %%solve the system

  sol= pdepe(m,@pdefun,@pdeic,@pdebc,x,t);

  % Extract the first solution component as u.
  u = sol(:,:,1); 

  % A surface plot 
  surf(x./rOuter,t./tau_chic_limiting,u-Tinf) 
  title(horzcat('Numerical solution computed with ',num2str(npoints), ' mesh points.'))
  xlabel('Nondimensional radial distance r/r_O_u_t_e_r')
  ylabel('Limiting Fourier number Fo')
  zlabel('Temperature rise T')
  set(gca, 'YScale','log') %if logspace used

  % A limiting solution profile.
  figure
  plot(x,u(end,:))
  hold on
  plot([r2 r2],[0 max(u(end,:))],'--')
  plot(x,0,'.')
  plot([r3 r3],[0 max(u(end,:))],'--')
  hold off
  title('Final timestep (Fo \rightarrow \infty) solution')
  xlabel('Distance r')
  ylabel('T_\infty')

  % --------------------------------------------------------------
  function [c,f,s] = pdefun(x,t,u,DuDx) %computes cfs

  %geometry and material info passed through nesting

  %material distribution of alpha, k

  if (x>=r2 && x<=r3) 
      alpha1=alpha_cond; 
      k=k_cond;
  else
      alpha1=alpha_ins;
      k=k_ins;
  end

  %define for pdepe
  c = 1/alpha1;

  %definition of flux term
  f = DuDx;

  %heating distribution (currently constant with t)

  rho_cond_max=rho_conductor; %for invariable resistivity
  sigma_cond=1/rho_cond_max; %definition of conductivity = 1/resistivity
  qvdot=ACfactor*(Jlimit^2)/sigma_cond;  %volumetric heating rate (of conductor) 

  %source terms: heating only in conductor 
  s_cond=qvdot/k_cond;
  s_ins=0/k_ins; 

  if (x>=r2 && x<=r3) 
      s=s_cond;
  else
      s=s_ins;
  end

  end

  % --------------------------------------------------------------
  function u0 = pdeic(x) %ICs
  u0 = 0.*x; %starts at zero relative temp
  end
  % --------------------------------------------------------------
  function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t) %BCs
  %exterior BCs only

  %constants passed through nesting

  %nb f is set as du/dx

  % %insulated LHS (inner)
  pl = 0;
  ql = -1*k_ins;

  %convective RHS (outer)
  pr = -1*h2*(ur-Tinf);
  qr = -1*k_ins; 

  end

  %end of nesting

  end