# How can I handle interface thermal boundary conditions?

I am trying to solve a system of PDE´s coupled by interface boundary conditions, this is my code:

g = 1;
k = 1;
f = 0.5;
Pe = 1;
pde =
{D[u[t, x, y], x, x] + k^2 D[u[t, x, y], y, y] ==
g D[u[t, x, y], t] + Pe (1 - (y/f)^2) D[u[t, x, y], x],
D[v[t, x, y], x, x] + k^2 D[v[t, x, y], y, y] == D[v[t, x, y], t]};

s = 1;
f = 0.5;
bc =
{u[t, 0, y] == Exp[-10000 t],
(u^(0,1,0))[t, 1, y] == Exp[-10000 t],
(u^(0,0,1))[t, x, 0] == Exp[-10000 t],
(u^(0,0,1))[t, x, f] == s (v^(0,0,1))[t, x, f],
u[t, x, f] == v[t, x,f],
u[0, x, y] == Exp[-10000 t],
(v^(0,1,0))[t, 0, y] == Exp[-10000 t],
(v^(0,1,0))[t, 1, y] == Exp[-10000 t],
v[0, x, y] == 1};

sol =
NDSolve[{pde, bc}, {u, v}, {t, 0, 6}, {x, 0, 1}, {y, 0, 1},
Method ->
{"MethodOfLines",
"SpatialDiscretization" ->
{"TensorProductGrid", "DifferenceOrder" -> "Pseudospectral"}}]


But the following error message appears:

NDSolve::bcedge: Boundary condition (u^(0,0,1))[t,x,0.5]==(v^(0,0,1))[t,x,0.5] is not specified on a single edge of the boundary of the computational domain. >>

Is here any way to fix this?

• Welcome to Mathematica.SE! I suggest the following: 0) Browse the common pitfalls question 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Read the faq! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! Oct 16, 2015 at 6:12
• You have things like v^(0,0,1) and u^(0,1,0) - which aren't valid Mathematica code. Oct 16, 2015 at 6:46
• @JasonB, when derivatives are copied as plain text they are represented exactly like so. It's likely valid syntax in the notebook. Oct 16, 2015 at 15:04

This is more an extended comment than an answer. If I have not made a trivial mistake, the boundary conditions in executable form are

bc =
{   u[0, x, y] == 0,
v[0, x, y] == 0,

u[t, 0, y] == Exp[-10000 t],
(D[v[t, x, y], x] /. x -> 0) == Exp[-10000 t],

(D[u[t, x, y], x] /. x -> 1) == Exp[-10000 t],
(D[v[t, x, y], x] /. x -> 1) == Exp[-10000 t],

(D[u[t, x, y], y] /. y -> 0) == Exp[-10000 t],

u[t, x, f] == v[t, x, f],
(D[u[t, x, y], y] /. y -> f) == s (D[v[t, x, y], y] /. y -> f)};


Running the code now produces the error described in the question. One course of action is to replace {y, 0, 1} by {y, 0, f} in NDSolve, with the intent of solving the PDE over the reduced range and then using that solution as boundary conditions to solve the PDE for {y, f, 1}.

However, there are other problems. First, the boundary condition for v at y -> 0 is missing and must be provided. (I reordered the boundary conditions to make this more apparent.) Second, several of the boundary conditions are inconsistent at boundary intersections. For instance, the first boundary condition in bc at x -> 0 is

u[0, 0, y] == 0


while the third boundary condition at t -> 0 is

u[0, 0, y] == Exp[0]


While such inconsistencies typically are not fatal, they can cause inaccuracies and should be fixed too. Finally,

"DifferenceOrder" -> "Pseudospectral"


is more appropriate for periodic grids and slows the computation somewhat. I recommend that the entire Method option be deleted. After these issues are resolved, it may be possible to tackle the most difficult problem that I have observed, namely that NDSolve does not seem to like the shared boundary conditions.

D[u[t, x, y], y] /. y -> f) == s (D[v[t, x, y], y] /. y -> f), u[t, x, f] == v[t, x, f]


which cause the error message,

NDSolve::bcsol: Could not solve for equations at boundary points from the boundary conditions. >>


I believe that this issue too can be overcome, probably by splicing u and v into a single function. However, it is pointless to try until the issues described above are resolved.