I need to know how to put the s variable of pdefun (from pdepe solver) as a function of x and time

Question

JOAN PERE PONSETI 2017-1-6

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/319350-i-need-to-know-how-to-put-the-s-variable-of-pdefun-from-pdepe-solver-as-a-function-of-x-and-time

评论： JClarcq 2018-2-19

Hi, I need help to solve a proble I'm in for a long ago. The thing is that I can't find any doc explaining how to set the variables [c,f,s] of pdefun (in particular the 's') as a function of the point of the 'xmesh'. It is to solve the heat transfer and difussion equation that is rho*cp*du/dt-k*d/dx(du/dx)=g.

from mathworks blogs and answers, I know that c=rho*cp (which is constant with time and x) f=k (also constant), but g is power generation, and that is not constant with time neither with x, the power is generated the first time values and just in a given set of points of xmesh.

Here I paste a part of the code where I state the pdefun:

function[c,f,s] = pdef(x, u, t, DuDx)
      c = rho*cp;
      f = k*DuDx;
      i=1;
      while(i<=length(x))
          if(x(i)<=percent*radi)
              s = (1/(h*pi*radi)^2)).*interp1(t40,Pot40,t,'linear','extrap');
          else
              s = 0;
          end
          i=i+1;
      end
end

Please, if anyone could tell me if this is well developed. The s variable I mean.

8 个评论
显示 6更早的评论隐藏 6更早的评论

JOAN PERE PONSETI 2017-1-11

Ok, thanks a lot. I will try.

JOAN PERE PONSETI 2017-1-11

在 MATLAB Online 中打开

Hi again, I implemented this code to see if now it works but it gives me this error message back.

Not enough input arguments.

Error in CLDissipation1D/pdef (line 375) if(x==xmesh(1))

Error in pdepe (line 246) [c,f,s] = feval(pde,xi(1),t(1),U,Ux,varargin{:});

Error in CLDissipation1D (line 340) sol = pdepe(m, @pdef, @pdeic, @pdebc, xmesh, t40);

function[c,f,s] = pdef(x,t,DuDx,S,xmesh,tspan)
c=rho*cp;
f=k*DuDx;
%  find out the volume we need the previous point too
i=1;
encontrado=0;
while(encontrado==0)
    if(x==xmesh(1))
        vol=inf;
        encontrado=1;
    elseif(x==xmesh(i))
        vol=(pi*(x^2-xmesh(i-1)^2)*h);
        encontrado=1;
    end
    i=i+1;
end
% repeat the same loop to find the power at time t
i=1;
encontrado=0;
while(encontrado==0)
      if(t==tspan(i))
          pot=Pot40(i);
      end
      i=i+1;
  end
% Finally we define the variable s as the ratio of power/volum
s=pot/(vol*pin);
% ========================================================================= % defines the Initial Conditions
end

Thanks a lot again!

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2017-1-11

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/319350-i-need-to-know-how-to-put-the-s-variable-of-pdefun-from-pdepe-solver-as-a-function-of-x-and-time#answer_250213

1. pdepe calls pdef for arbitrary times, not only those specified in tspan.

2. sol = pdepe(m, @(x,t,DuDx)pdef(x,t,DuDx,S,xmesh,tspan), @pdeic, @pdebc);

3. Comparisons for equality like if x==xmesh(1) probably won't work because of rounding errors.

Best wishes

Torsten.

9 个评论
显示 7更早的评论隐藏 7更早的评论

JOAN PERE PONSETI 2017-1-13

在 MATLAB Online 中打开

And now, this. Here the problem is in pdeic??

Error using CLDissipation1D/pdeic Too many input arguments.

Error in pdepe (line 229) temp = feval(ic,xmesh(1),varargin{:});

Error in CLDissipation1D (line 188) sol = pdepe(m, @(vol)pdef,@pdeic, @pdebc, xmesh, t40,[], vol);

sol = pdepe(m, @(vol)pdef,@pdeic, @pdebc, xmesh, t40,[], vol);
u=sol(:,:,1);
%========================================================================= % defines the properties of the PDE function
function[c,f,s] = pdef(x,t,u,DuDx,vol)
c=rho*cp;
f=k*DuDx;
s=interp1(t40,Pot40,t,'linear','extrap')/(pin.*interp1(xmesh,vol,x,'spline','extrap'));
% ========================================================================= % defines the Initial Conditions
end
function[u0] = pdeic(x)
u0 = tinitial;
end % ========================================================================= % defines the Boundary Conditions
function[pl,ql,pr,qr] = pdebc(xl, ul, xr, ur, t)
%for cylindrical pl and ql are 0, because of symmetry. No flux is present
%at the center of the sphere.
pl = 0;
ql = 0;
pr = ur-interp1(t40,Temp40,t,'linear','extrap');
qr = 0;
end

JOAN PERE PONSETI 2017-1-21

在 MATLAB Online 中打开

Thanks a lot! I've solved this by stating the variables as global. But now, I have another problem with the solution I get. The thing is that, there is a power dissipation at the center, and I have the initial conditions and the boundary ones. Logic says that if there is power dissipation inside, the temperature within the core must be always higher than the one outside (at the same time for both temperatures). But what I get is the graph below, and furthermore, the blue-crossed line of the graph doesn't start at the 300 K (aprox) which is the initial condition.

I give you the code I have used and let's see if you can figure out what's the error.

Thanks!

Joan Pere

It looks like the differential is not working properly.. don't know ehat happens

    t40=times{var};
        Temp40=temperatures{var};
        Pot40=potencies{var};
    % Assignem els vectors a variables que es mostraran en el Workspace
    assignin('base','t40',t40);
    assignin('base','Temp40',Temp40);
    assignin('base','Pot40',Pot40)
    % figure
    % plot(t40,Pot40);
    % 
    % figure
    % plot(t40,Temp40);
    % capacitor properties
    k = 1.5; % W/(m K) 30 Alumina %aluminum 237 
    rho = 3000; % alumina density 4000 Aluminum 2300 
    cp = 900; % J/K/kg Aluminum 910 Alumina 880
    %capacitor geometry
    h=0.01;
    radi=0.005;
    percent=0.25;
    tinitial = Temp40(1); % intial temp throughout the tile
    m = 1; % assume cylindrical symmetry
    xmesh = linspace(0, radi, 200); % 100 punts a discretitzar
    assignin('base','xmesh',xmesh)
    global vol
    for i=1:length(xmesh)
         if(i==1)
             vol(i)=pi*h*(radi)^2;
         else
             vol(i)=inf;
         end
     end
    assignin('base','vol',vol)
    sol = pdepe(m, @pdef,@pdeic, @pdebc, xmesh, t40);
    T1=sol(:,:,1);
    assignin('base','T1',T1)
    % A surface plot is often a good way to study a solution.
    figure;
    subplot(2,1,1)
    surf(xmesh,t40,T1);
    title(sprintf('SOLUCIÓ NUMÈRICA AMB 200 PUNTS.\nArxiu %s.', VCC40));
    xlabel('Distància al centre (m)');
    ylabel('Temps (segons)');
    % A solution profile can also be illuminating.
    subplot(2,1,2)
    plot(t40,T1(:,1),'x');
    hold on
    plot(t40,T1(:,end),'+');
    hold on
    plot(t40,Temp40,'o');
    legend('Temperatura al centre','Temperatura externa','Condicions Incials','Location','best');
    title(sprintf('PERFIL DE TEMPERATURES EN 2 PUNTS (x=0 i x=R).\nArxiu %s.', VCC40));
    xlabel('Temps (segons)');
    ylabel('Temperatura (Kelvin)');
    %========================================================================== % defines the properties of the PDE function
    function[c,f,s] = pdef(x,u,t,DuDx)
    c=rho*cp;
    f=k*DuDx;
    s=(1/interp1(xmesh,vol,x,'linear','extrap')).*interp1(t40,Pot40,t,'linear','extrap');
    end
    % ========================================================================= % defines the Initial Conditions
    function[u0] = pdeic(x)
    u0 = tinitial;
    end
    % ========================================================================= % defines the Boundary Conditions
    function[pl,ql,pr,qr] = pdebc(xl, ul, xr, ur, t)
    %for cylindrical pl and ql are 0, because of symmetry. No flux is present at the center of the sphere.
    pl = 0;
    ql = 0;
    pr = ur-interp1(t40,Temp40,t,'linear','extrap');
    qr = 0;
    end

The volume is needed because my source term is the power density (Power/Volume) then the power is a different one each time (that's why I have to interpolate) and the volume I state that the volume dissipation is a quarter of the radius (percent=0.25)

JOAN PERE PONSETI 2017-2-24

在 MATLAB Online 中打开

Another thing, when I write the next code % Specify Coefficients

f = @(region,state) interp1(t40,Pot40,state.time,'linear','extrap');
specifyCoefficients(pdem,'m',0,'d',rho*cp,'c',k,'a',0,'f', f);

and I set a breakpoint at f, it always says that state structure is 0 at all its elements but it doesn't return any error until the fifth or sixth step... don't know why...

thanks!

JClarcq 2018-2-19