Using pdepe to solve heat diffusion between two laysers

To learn to use pdepe to solve the heat diffusion between two layers, I started with using pdepe to solve the classic problem "two semi-infinite bodies in contact," whose theoretical solution predicts that the interface temperature would reach to a constant temperature at the moment of contact and remain constant throughout the contact period. However, the interface temperature calculated by pdepe is time dependent though the value is quite close to the theoretical one. I appreciate if anyone could advice me why pdepe doesn't give the "correct answer"? Is it because there is sth wrong with my code?

The heat equation that I'm trying to solve is: rho*c*du/dt=k*d2u/dx2

and my code is as follows:

==========================

function pdex4

m = 0;

x=[-100, -0.2, -0.1, -0.001,-0.0001,-0.000001, 0, 0.000001,0.0001,0.001, 0.1, 0.2, 100];

% -100<x<0, object 1 % 0<=x<100, object 2

t = [0 0.001 0.01 0.1 0.5 1 1.5 20];

sol = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t);

u1 = sol(:,:,1)

% --------------------------------------------------------------

function [c,f,s] = pdex4pde(x,t,u,DuDx) %governing equation x_interface=0;

if x < x_interface %object 1

    c=rho1*c1;
    f = k1*DuDx; 
    s = 0;

else %object 2

    c=rhos*cs;
    f = ks*DuDx; 
    s = 0;

end

% --------------------------------------------------------------

function u0 = pdex4ic(x)

x_interface=0;

if x < x_interface %object 1 initial condition

    u0=u1i;

else

    u0=u2i;           %object 2 initial condition

end

% --------------------------------------------------------------

function [pl,ql,pr,qr] = pdex4bc(xl,ul,xr,ur,t)

% BC for semi-infinite bodies: Same as the initial conditions at two ends.

pl = ul(1)-u1i; ql = 0; pr = ur(1)-u2i; qr = 0;

==========================

Answer 1

[Torsten]

You set the initial temperature at the interface to u2i - so you can't expect that pdepe knows better.

Best wishes

Torsten.

[asphalt]

Dear Torsten,

Thank you very much for your reply. Actually I'm quite confused at how to set the conditions (eg. c,f,s and u0) at the interface. Since there is only one Xmesh point for the interface(x=0 in my case), should the interface have the properties of object 1 or object 2 or none of them? Is skipping the interface (x=0) when setting c,f,s and initial conditions (see below) a more proper way when dealing with the interface?

==

if x < x_interface %object 1 initial condition

u0=u1i;

elseif x > x_interface

u0=u2i; %object 2 initial condition

end

==

Thank you, asphalt

[Torsten]

It's not possible to skip the conditions for the interface.

Usually, the properties are set by the arithmetic mean of the properties from the right and from the left.

By the way: Could you explicitly state the problem you are trying to solve and for which you think pdepe does not give the correct answer ?

Best wishes

Torsten.

[asphalt]

Dear Torsten,

Thank you very much for the response again.

I'm trying to solve the temperature responses at the interface when two objects are brought into sudden contact.

In my code I assume that the contact is perfect so that the interface temperature at both side is the same and I set the thickness of the objects to be very large (x=100). Under this situation, the theoretical solution indicates that the interface temperature should change to an intermediate temperature of the initial temperatures of objects 1&2 (e.g. 24.5) at the moment of contact and remain constant throughout the contact. However, the solution given by pedepe will first reach a lower temperature (e.g. 24) and then increase to a higher one (e.g. 24.8) and then decrease again to the lower temperature (24.2) as time increases. Since this solution is time dependent, I'm not sure whether it is because there is some problem in the way that I set the problem or it is a property of PDEPE.

The other question that I want to ask is whether it is suitable to use a single PDE to model two objects brought into sudden contact. In my real problem, there will be a thermal contact resistance at the interface and the interface temperature at the two sides will be different. In this case, a single PDE seems not be able to solve this problem because I can't apply BC (flux=(T1-T1)/R)at the interface. Would using 2 PDEs (one for each object) a better choice? Or should I use PDE toolbox for this problem?

Best, asphalt

[Torsten]

It will be better to solve both problems in two separate domains.

The problem you tried above will have transmission conditions

T1 = T2
lambda1*dT1/dn = lambda2*dT2/dn

at the interface during the calculation.

pdepe is not suited for this purpose. I don't know if problems with contact resistance can explicitly be treated with the PDE Toolbox. I'd discretize the equations in space (1d), incorporate the transmission condition and solve the resulting system of ordinary differential equations using ODE15S, e.g.

Best wishes

Torsten.