How to use the PDE Toolbox combined with the script?

Hello, is there anybody who can help me with my question? I need to solve a 2D heat diffusion equation PDE, I have the boundary conditions and initial conditions set in an excel file, I also have them in a script. So, as the boundary conditions vary with time but I haven't got any expression for that but the file, I need to import this data from the script or excel file to PDE toolbox.

The thing is the following:

- Boundary Conditions: ExternalTemperature(t) is a vector of a given length containing the external temperatures of a cylinder, so at each time, the bc's change.

- Initial Conditions: IC(0), is the first temperature in that previous vector. ExternalTemperature(1).

There are two zones at the cylinder, and outer one (external crown) which doesn't generates power, so the source term is 0, and an inner one which does generate power, and has a given value for the source term. So, I need to set the coefficients as a function of time (power is a function of time) and space (x-y) since the points outside the inner cylinder doesn't generate power.

Thanks in advance,

Joan Pere!

PS: this is for a final degree project, so any help is gratefully welcomed!

[Sam McDonald]

Joan Pere,

It is difficult to address your question without additional information. Could you provide some more detail on where you are getting stuck? Is it loading the data into MATLAB, or using it in the PDE Toolbox to define the boundary conditions? Do you have some example code you are writing or modifying to solve this problem?

In the meantime, if you have not already, take a look at these links from the MathWorks documentation:

https://www.mathworks.com/help/pde/ug/steps-to-specify-a-boundary-conditions-object.html https://www.mathworks.com/help/pde/ug/boundary-object-with-nonconstant-conditions.html

[JOAN PERE PONSETI]

Hi, thanks for the response. I have decided to go back and use the matlab functions like solvepde for 2D resolution. But now, I have a problem when it comes to set the geometry. I need a layer geometry but matlab returns an error when running the code. The error is the following:

Error using CoreDissipation2D (line 269)
Invalid geometry detected. Edges overlap or intersect at non-end points.

The code is simple: I only need concentrical circles, and don't know why it is seen that Matlab cannot execute it.

ncapes=6;
R2=[1;0;0;radius]; % contorn exterior (transferència de temperatura)
% Per a crear les geometries de les altres capes, necessitarem ncapes-1
% circunferències més, per això definrem una matriu per establir les
% circunferències internes.
% nint=ncapes-1;
I1=[1;0;0;radius/ncapes]; % primer contorn interior
I2=[1;0;0;2*radius/ncapes]; % segon contorn interior
I3=[1;0;0;3*radius/ncapes]; % tercer contorn interior
I4=[1;0;0;4*radius/ncapes]; % quart contron interior
I5=[1;0;0;5*radius/ncapes]; % cinquè contron interior
gd2 = [R2,I1,I2,I3,I4,I5];
sf2 = 'R2+I1+I2+I3+I4+I5';
ns2 = char('R2','I1','I2','I3','I4','I5')';
g2 = decsg(gd2,sf2,ns2);

What I understand from the error is that the circles overlap or intersect each other, but don't really think this is the problem, because I have already tried another geometry similar to this but with only two circles.

Thanks a lot,

Joan Pere.

Answer 1

[Alan Weiss]

It is possible that you are getting bit by a bug. Sorry. I believe that you will have no problems if you create your circles using a geometry function. For example,

function [x,y] = multicircle(bs,s)
% Create circles centered at (0,0) using four segments per circle

switch nargin
    case 0
        x = 6*4; % four edge segments per radius
        return
    case 1
        A1 = repmat([0,pi/2,pi,3*pi/2],1,6); % start parameter values
        A2 = repmat([pi/2,pi,3*pi/2,2*pi],1,6); % end parameter values
        A3 = repmat(1:6,4,1);
        A3 = A3(:)'; % region label to left
        A4 = repmat([2:6,0],4,1);
        A4 = A4(:)'; % region label to right
        A = [A1;A2;A3;A4];
        x = A(:,bs); % return requested columns
       return
    case 2
        if numel(bs) == 1
            bs = bs*ones(size(s));
        end
        ind = floor((bs-1)/4) + 1; % circle number from 1 through 6
        x = ind.*cos(s);
        y = ind.*sin(s);
end

I tested this geometry using the following script:

model = createpde();
geometryFromEdges(model,@multicircle);
generateMesh(model);
pdegplot(model)
pdemesh(model)
axis equal

I hope that this helps. Sorry about the bug.

Alan Weiss

MATLAB mathematical toolbox documentation

[JOAN PERE PONSETI]

Ok, I will try this one. Am I supposed to give the function the variables bs and s. So, what is the format of these variables? I guess it is the radius of the circles, isn't it?

Thanks,

Joan Pere

[Alan Weiss]

Don't worry about bs and s, MATLAB will pass those variables to your geometry function. I suggest that you try the function exactly as I gave it and see how well it works for you.

The explanation of the syntax is in the link I gave in the first paragraph of my previous answer.

Alan Weiss

MATLAB mathematical toolbox documentation

[JOAN PERE PONSETI]

Hello, I'm facing now another issue. The thing is that in specifyCoefficients function I need to set f coefficient as a function of x and y coordinates and time. I'm using the next function:

function [f]=f1(region,state) % =========================================== Defines the f coefficient as a function of the space (region.x ; region.y) and time (state.time)

k=1;
while(k<=length(region.x)) 
% length(region.x) should be the number of nodes, isnt it?
% the use of k index is to define f at each node, but I need to take into account the time too

       if(region.subdomain==1)
           f(1,k)=interp1(t40,Pot40,state.time,'linear','extrap')/(volumeD.*h);
       else
           f(1,k)=0;
       end

   k=k+1;

end

end

It is supposed to compute the f coefficient depending on x and y, and time too. From Matlab documentation, region.x region.y region.z region.subdomain are row vectors, but to make sure I define a variable a (for example) and give it the value of region.x. What I have in return is a single double value (0 always). So, is it a vector, a single value...? The function is supposed to return f as a vector, isn't it? But, as I need the f coefficient in space and time, wouldn't it be kind of matrix, for each time, and node I will have a given value of f coefficient.

Thanks a lot,

I'm in a hurry with this because I am supposed to compute this heat transfer equation for a two layered and multi layered geometry. And none of them gives me results back, as I expect them to be.

Joan Pere

[Alan Weiss]

I'm not sure that you have a firm idea about what a coefficient function is. MATLAB passes the region and state variables, and your function returns a vector of values corresponding the coefficient at each of the points in the region structure, possibly taking into account the state structure. (In general, your function returns a matrix, but because you have a scalar f, it returns a row vector.) See the documentation.

I do not understand your code because I do not know what t40, Pot40, volume, or h represent. You also do not seem to be using the region.x or region.y fields, or the region.subdomain field. You almost certainly want to use those fields if your f coefficient is not a constant.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation

[JOAN PERE PONSETI]

Hi, I've been working on this but don't get any results. The thing is that, I haven't got a continuous function that defines the coefficients like f=y^2*sin(t), in this case, it'd be direct to state :

f = @(region,state) region.y.^2*(sin(state.time)).

But I have a non-continuous function, that is, I haven't got a mathematical expression for that, so I have to use a while, and depending on the x, y and t variables (region.x, region.y and state.time) I have to set the coefficients. Have you any idea of how to solve it??

Thanks a lot!

[Alan Weiss]

Of course, I don't know what your function is. But if you can write a MATLAB function that gives you the coefficient value depending on x, y, and t, then you can use that function as the basis of your coefficient function. And if you cannot write a MATLAB function that gives you the coefficient value, then I don't see how you can expect to define the problem in MATLAB.

For a more useful response, give more details about your coefficient function, and I can then try to show how to include it in PDE Toolbox format.

Alan Weiss

MATLAB mathematical toolbox documentation

[JOAN PERE PONSETI]

Ok, I have a cilynder (capacitor) which is modelled as a layered section. This is, n layers and each layer is formed with 3 sublayers (Aluminum-Dielectric-Aluminum). It is repeated n times, as much as layers I have. To sum it up, the odd layers (aluminum) must have the aluminum coefficients and the even ones must have the dielectric material coefficients. There is an external layer, an insulating cover which doesn't generate power, so f would be 0.

The inner layers, must have the following coefficients:

aluminum layers: m=0; d=rho*cp; c=k; a=0; f=power/volume;

dielectric layers: m=0; d=rhodielectric*cpdielectric; c=kdielectric; a=0; f=power/volume;

insulating layer: m=0; d=rhoinsulating*cpinsulating; c=kinsulating; a=0; f=0;

where all rho, cp and k are scalars, thermal properties of the materials.

I am wondering if it is possible that, with one function, I can have three outputs, so I can set d c and f using the same function, like.

function [d,c,f]=coefficientsfunction(region,state)

Is the idea clear? Don't know if I explained myself the best.

Thanks a lot!

[Alan Weiss]

Thank you for providing some details. I am going to assume that you have a 2-D model with subdomains representing your various "layers." If you have a 3-D model then I am afraid that you are out of luck, as, currently, there are no subdomains for 3-D.

For your f coefficient, you can try something like the following. I am assuming that the insulating layer is region #1, so modify this script as necessary if I am wrong about your geometry, and I am assuming that the even layers are dielectric. I also assume that you have a constant pv = power/volume that you can pass to the function.

function f = fcoef(region,state,pv)
L = length(region.x);
f = zeros(1,L);
odds = mod(region.subdomain,2) == 1;
evens = mod(region.subdomain,2) == 0;
f = pv*odds + pv*evens;
oneregion = region.subdomain == 1;
f(oneregion) = 0;

In your coefficient setting statement, pass @(region,state)fcoef(region,state,pv).

Similarly, for the d coefficient:

function d = dcoeff(region,state,rho,cp,...
  rhodielectric,cpdielectric,rhoinsulating,cpinsulating)
L = length(region.x);
d = zeros(1,L);
odds = mod(region.subdomain,2) == 1;
evens = mod(region.subdomain,2) == 0;
d(evens) = rhodielectric*cpdielectric;
d(odds) = rho*cp;
oneregion = region.subdomain == 1;
d(oneregion) = rhoinsulating*cpinsulating;

Pass this in your coefficient setting statement as @(region,state)dcoeff(region,state,rho,cp,rhodielectric,cpdielectric,rhoinsulating,cpinsulating)

You can write a c coefficient along the same lines.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation

[JOAN PERE PONSETI]

Don't know exactly if I can define the functions like this, I have no subdomains defined. I have only a circle with radius r and then, depending on the distance (not comparing with subdomains) I set ones or others coefficients.

I have tried with this one:

function f=fcoeff(region,state)

      l=length(region.x);
      f=zeros(1,l);
      a=hypot(region.x,region.y); % distance to the centre of all nodes

      for i=1:l

          % Let's check in which layer these distances are.
          j=1;
          found=0;
          while(j<length(vectorlimits) && found==0)

              if(a(l)>=vectorlimits(j) && a(l)<=vectorlimits(j+1))

                  layer=j; % remember that j=1 means the most inner layer (Aluminum)
                  layervolume=pi*(vectorlimits(j+1)^2-vectorlimits(j)^)*h;

                  % if the ayer is the last one, then, it is the insulating one
                  if(j==length(vectorlimits)-1)
                      f(1,l)=0; % There is no power dissipation
                  else

                      % Once that we know it is not the last layer, we must
                      % check if it is aluminum (odd layers) or dielectric
                      % (even layers)

                      if(mod(layer,2)~=0) %odd layer, aluminum properties, we distribute the power percentualy depending 
                          % on the volume of the layer, the higher the volume, the hiher the power (equal power density)
                          f(1,l)=interp1(t40,Pot40,state.time,'linear','extrap').*volumspercent(layer)/layervolume;
                      else % even layer
                          f(1,l)=interp1(t40,Pot40,state.time,'linear','extrap').*volumspercent(layer)/layervolume;
                      end
                  end

              else
                  j=j+1;
              end
          end
      end

end

I know that the terms where f is not 0, can be combined with just one line since the material, for the source term, is irrelevant, but I'd rather have like this to implement the same function, if I can, to the other coefficients.

The line where I set a=hypot(...,...); Theoretically it returns a vector with the distance to the center isn't it? I tell you because in past tries I have been given by matlab that 'a' or region.x or region.y are 0.

Thanks a lot! Joan Pere.

[Alan Weiss]

I think that you are making a serious mistake in not defining subdomains. The mesh will not respect the boundaries that are in your head, and therefore the elements can overlap the regions. I mean your mesh will include points from more than one type of region, no matter how fine your mesh.

Also, I have trouble understanding your code because it seems to be written with the assumption that the regions are in some kind of order. I don't see how you are passing in the vectorlimits data, or t40 or Pot40. I think that this line is erroneous:

if(a(l)>=vectorlimits(j) && a(l)<=vectorlimits(j+1))

You should decide whether you want to write code in a vectorized fashion or in a loop. It is faster to write vectorized code, but until you really know what you are doing, write loops. I assume that vectorlimits is monotone increasing.

for i = 1:L % I prefer upper-case L, l looks like 1
  indx = find(vectorlimits - a(i) > 0,1) % find the band
  if mod(indx,2) = 1 % do your calculation here
    f(i) = ...
  else
    f(i) = ...
  end
end

I hope that this helps,

Alan Weiss

MATLAB mathematical toolbox documentation

[Robin]

Hello Alan,

I am having a similar problem. In my case i want to create several circular structures in different locations. My problem is that I don't know how to modify your code. If I'm understanding the code correctly, I just wanted to adapt the last part under case 2. So my idea was to shift the cos or sin function in my desired direction.

x = ind.*cos(s + 1);

y = ind.*sin(s - 1);

But this results in the following error: Error using test_coefficient_def_for_geometry

Meshing failed due to invalid geometry. Each face must have a unique face ID.

Could you elaborate on this?