How can I simulate 2 materials in pde toolbox script?

In order to simulate 2 materials in a heat diffusion equation I define 2 different geometries, but then how I define the different coefficients for each material? One has heat source and the other hasn't. I am doing it as shown bellow, but the results are not what I would expect, so I am wondering if I am doing it right. There is another way to do that? Do I need to define 2 diff. equations instead?

Thank you very much in advance,

Here is my script:

numberOfPDE = 1;
model = createpde(numberOfPDE);
r1 = [3,4,0,0,0.5*10^-6,0.5*10^-6,-160*10^-9,-1.16*10^-6,-1.16*10^-6,-160*10^-9];
r2 = [3,4,0,0,0.5*10^-6,0.5*10^-6,0,-160*10^-9,-160*10^-9,0];
gdm = [r1; r2]';
g = decsg(gdm,'R1+R2',['R1'; 'R2']');
geometryFromEdges(model,g);
specifyCoefficients(model,'m',0,'d',5,'c',16,'a',0,'f',0,'face',1);
specifyCoefficients(model,'m',0,'d',5,'c',50,'a',0,'f',2*10^18,'face',2);
applyBoundaryCondition(model,'edge',1,'r',298);
applyBoundaryCondition(model,'edge',2,'q',2000,'g',596000);
applyBoundaryCondition(model,'edge',3,'g',0);
applyBoundaryCondition(model,'edge',4,'g',0);
applyBoundaryCondition(model,'edge',6,'g',0);
applyBoundaryCondition(model,'edge',7,'g',0);
hmax = 2*10^-8; % element size
msh=generateMesh(model,'Hmax',hmax);
tinitial = 0;
tfinal = 2*10^-8;
interval = 100;
tlist = linspace(tinitial,tfinal,interval);
u0 = 298; % initial temperature in the geometry 
setInitialConditions(model,u0);
results = solvepde(model,tlist);
u = result.NodalSolution;
%Plot the solution at time t = tfinal.
figure;
pdeplot(model,'xydata',results.NodalSolution,'contour','on')

[michio]

It seems to me that you are doing alright in specifying different coefficients for different sub-domains. I assume two different materials are represented by the two sub-domains.

specifyCoefficients(model,'m',0,'d',5,'c',16,'a',0,'f',0,'face',1);
specifyCoefficients(model,'m',0,'d',5,'c',50,'a',0,'f',2*10^18,'face',2);

It may help us help you if you describe more about your issues, your problem settings, what part of the result is not what you expect.

[Laura Monreal]

Hi Michio,

I changed my name because I had problems with the other sesion. Thank you for you answer. So, in fact my f coefficient for the subdomain 2 is a cosinus that I do not see in the simulation, but i see it when I simulatet the "material" 2 byitself. I think I am not getting the result I spect because of the mesh, there is a way to introduce different meshing for different subdomains? I saw that you can not implement refinemesh usind a PDEmodel, so how should I do it?

Thank you in advance

Answer 1

[michio]

Thanks for your clarification. As you have already made use of it, but sub-domain support is available for 2D problems, so yes you can simulate 2 materials using PDE Toolbox.

If your issue is really caused by meshing, you can try legacy workflow that supports refinement. This entry might be of your help. How can I refine a subdomain in the PDE Toolbox mesh generation tool?

I didn't understand the details since you mentioned "my f coefficient for the subdomain 2 is a cosinus". If your f coefficient needs to be expressed as some function of space, you can try defining f coefficient in a functional form. The related doc page is f Coefficient for specifyCoefficients

[Laura Pare]

Hi Michio,

Thank you for your help, I am already defining coefficients in a functional form, I just simplify the code I copied in order to see it easier. About the refine subdomains: I already tried the link you are sending me and it works, but the problem is that then when I try to solve the pde it tells me that the geometry doesn't have a mesh...and I have read that you can't refine mesh in a PDEModel: https://es.mathworks.com/matlabcentral/answers/236357-refine-mesh-in-pdemodel

So, I am not sure if I am doing something wrong or you can just not refine a subdomain and then solve the equation...

Thank you

[michio]

Hi Laura,

You are correct, subdomain refinement is not currently available for PDEModel, the recommended workflow, but I believe it's doable in the legacy workflow.

The URL I put in the above answer was wrong. I apologize. Please refer to https://jp.mathworks.com/matlabcentral/answers/302250-how-can-i-refine-a-subdomain-in-the-pde-toolbox-mesh-generation-tool

An example of a legacy workflow is to use hyperbolic function (R2016a)

The page on hyperbolic is depreciated in R2016b, so I put a direct link to the one from R2016a. You might need to login to your Mathworks account to visit.

[Laura Pare]

Hi Michio,

Thank you very much for your answer. It was very helpful, I am doing it in the legacy workflow and even without refining the mesh differently for the different subdomains I see results that make much more sense to me. But the only thing I don't know how to do now is to express a constant coeffitient for different subdomains.

I have this:

   u = parabolic(u0,tlist,model,'ccoeff(sd)',a,f,'dcoeff(x,y,t,u,sd)');

And the d is function like that:

if true

function d = dcoeff(x,y,t,u,sd)
d = 3950*(92.4+3.8*10^(-2).*(u)-2.2*10^6.*(u).^(-2))/101.9612*1000; % d on subdomain 1 (znO)
r = '(sd == 2)'; % subdomain 2
d2 = 5600*(47.8+6.1*10.^(-3)*(u)-7.6*10^5.*(u).^(-2))/81.3794*1000; % coefficient on subdomain 2 (Alumina)
d(r) = d2(r); % d on subdomain 2
end

But the c coeffitient is constant in both subdomains but differwent between them. I write this, but it tell me Index exceeds matrix dimensions.In expression: ccoeff(sd) when evaluating pde coefficients.

if true
function c = ccoeff(sd)
c = 50; % d on subdomain 1 
r = '(sd == 2)'; % subdomain 2
c2 = 16; % coefficient on subdomain 2 
c(r) = c2(r); % d on subdomain 2
end

Do you know why?

Thank you very much

Laura

[michio]

It would be great if I could reproduce the issues in my end, but I am guessing c(r), indexing in to a constant, causes the problem. Why don't you try

function c = ccoeff(sd)
c = 50*ones(size(sd)); % d on subdomain 1 
r = '(sd == 2)'; % subdomain 2
c2 = 16*ones(size(sd)); % coefficient on subdomain 2 
c(r) = c2(r); % d on subdomain 2
end

instead. Let me know how it goes :)

[Laura Pare]

It works perfectly,

Thank you very much for your help Michio!

[michio]

You are welcome!

Answer 2

[Alan Weiss]

You might like the plot better if you plot the end time rather than the initial time:

 pdeplot(model,'xydata',results.NodalSolution(:,end))

I am not sure that you have scaled things correctly, because I see very little difference between this plot and the plot at the second time step:

 figure;pdeplot(model,'xydata',results.NodalSolution(:,2))

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation

[Laura Pare]

Hi Alan,

Thank you for your answer. You are right about the time, I copied that part of the script wrong.

But then, is that the way to simulate two different materials with different coefficients?

I also saw something separating the coefficients by ! , but I didn't get how to do it: http://es.mathworks.com/help/pde/ug/scalar-coefficients-in-string-form.html?searchHighlight=subdomains

Shoudn't I simulate 2 differential equations then?

Thanks,

Answer 3

[Moncef]

Hello,

please how to introduce the c coeffecient which is has matrix form:

c = [x-y, 0; 1, x.^2]

using "specifyCoefficients" matlab command?

With best regards