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

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.

[Nikhil Sreekumar]

Hi Joan,

I had some doubts could you please clarify?

• Are you looking into the equation1) rho * Cp * (du/dt) = k * (d^2u/dt^2) + s? and 2) deltaT/deltat = (s*L) / (h * A)?
• If we are rewriting the equation 2) like(Potf - Poti) / (tf - ti) = (s * L) / (h * A)then with the available value of Potf, Poti, tf, ti, h and A, you can find the value of s right?

Thanks

Nikhil

[JOAN PERE PONSETI]

Hi, yes, my doubt is regarding to the equation you call 1), but variable u in my case is temperature T. The thing is that, matlab says that s is a function of x,u,t,dudx and I just want it to put s(x,t) because the source term depends of the position within the radius (x) and the time t, with the power. My source term is the power density generation, which is P/volume and the power depends on time and the volume depends on x. this is why my source term depends on x and time.

thanks!

JP

[Torsten]

The formulation

"s is a function of x,u,t,dudx"

means that s may depend on all of these variables, but of course it is possible that it depends on x and t only.

Best wishes

Torsten.

[JOAN PERE PONSETI]

Yes, I know, but my problem is that I don't know exactly how the anonymous function of pdefun works, I mean, if the pdepe calls for the pdefun each time in 'tspan' and each x in 'xmesh', if so, pdefun returns 's' as a single value (double) but don't know how to give the function the vectors of tspan and xmesh and then return that single value.

Thanks!

[JOAN PERE PONSETI]

The thing is that there is a velue for s in each time and each x, then, for a given x point there will be never the same value of s because it also changes on time. So, if I have a tspan of 20 intervals, and and xmesh of 50 points, s should be a matrix (50x20) which will have in each cell the value of the source term at each time and each x.

[Torsten]

Why do you think you have to supply tspan and xmesh in "pdefun" ?

From the documentation of pdepe:

pdefun computes the terms c, f, and s (Equation 1-3). It has the form 

[c,f,s] = pdefun(x,t,u,dudx)

The input arguments are scalars x and t and vectors u and dudx that approximate the solution u and its partial derivative with respect to x, respectively. c, f, and s are column vectors. c stores the diagonal elements of the matrix c (Equation 1-3).

So your code should read:

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

Of course, the variables percent, radi, h, t40 and Pot40 must be visible in "pdef".

Best wishes

Torsten.

[JOAN PERE PONSETI]

Ok, thanks a lot. I will try.

[JOAN PERE PONSETI]

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]

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.

[JOAN PERE PONSETI]

Ok, I have tried this, but matlab returns this error.

Not enough input arguments.

Error in CLDissipation1D>@(x,t,DuDx,vol,xmesh)pdef(x,t,DuDx,vol,xmesh) (line 296) sol = pdepe(m, @(x,t,DuDx,vol,xmesh)pdef(x,t,DuDx,vol,xmesh), @pdeic, @pdebc, xmesh, t40);

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

Error in CLDissipation1D (line 296) sol = pdepe(m, @(x,t,DuDx,vol,xmesh)pdef(x,t,DuDx,vol,xmesh), @pdeic, @pdebc, xmesh, t40);

sol = pdepe(m, @(x,t,DuDx,vol,xmesh)pdef(x,t,DuDx,vol,xmesh), @pdeic, @pdebc, xmesh, t40);

function[c,f,s] = pdef(x,t,DuDx,vol,xmesh)

c=rho*cp;

f=k*DuDx;

s=interp1(t40,Pot40,t,'linear','extrap')/(pin.*interp1(xmesh,vol,x,'linear','extrap'));

% ========================================================================= % defines the Initial Conditions
end

Thanks a lot!

[JOAN PERE PONSETI]

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

[Torsten]

Take a look at this page on how to pass extra parameters by the anonymous function method:

https://de.mathworks.com/help/optim/ug/passing-extra-parameters.html

Best wishes

Torsten.

[JOAN PERE PONSETI]

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)

[Torsten]

1. In the boundary part (@pdebc), pl=ql=0 can't be correct.

2. I don't see where "percent" is used in the calculation of the source term.

3. I don't understand your variable "vol" because you assign "inf" to all the inner grid points.

Best wishes

Torsten.

[JOAN PERE PONSETI]

1. pl and ql must be zero because there is no te flux into that point, which is the center of the cilindre.

2. the percent variable let's say is irrelevant here.

3. I see now that the volume variable code is not complete, I might have deleted some code lines. But the thing is that, i need the inner points to have a given value (power/volume, both of them are scalars) but i need the outter ones to be 0, thats why I set de volume to inf, since it divides and makes the value of the source term 0 (power always have a positive value).

Thanks!

[JOAN PERE PONSETI]

Hello again, I need some help in another code, similar to this one but now using the solvepde function. When I have to assign the coefficients. I guess that when I write the following code:

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

Matlab returs an error that says the following:

Error using pde.CoefficientAssignment/checkCoefFcnHdlArgCounts (line 531)
Function handle specifying a coefficient must accept two
input arguments and return one output argument.

As if the function doesn't return any value, is it parmited to use the state and region variables in a function like interp1??

Thanks a lot!!

Joan Pere

[JOAN PERE PONSETI]

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]

Hi Joan,

I have the same problem than you where my heatsource is a term temperature dependent via an interp1 of a material properties.

Have you been able to solve this inputs argument issue? Thanks,