# How to solve a PDE with Robin Boundary Condition inside considered region?

I'm trying to solve the heat diffusion equation in cylindrical coordinates. The main problem is that I would like to include the Robin Boundary Condition inside considered region in order to simulate the interface between two materials. Particularly, I would like to have the value of the derrivatiove on the internal boundary proportional to the function difference on both sides of the boundary (here z=0). Of course, the function should not be continous for z=0.

Is it possible in Mathematica 10? Below is my code. Thanks in advance!

Needs["NDSolveFEM"];
Rmax = 5; Lmax = 0.2;
WholeRegion = ToBoundaryMesh[
"Coordinates" -> {{0, Lmax}, {0, 0}, {0, -Lmax}, {Rmax,
Lmax}, {Rmax, 0}, {Rmax, -Lmax}},
"BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {4, 5}, {5,
6}, {1, 4}, {2, 5}, {3, 6}}]}
];
WholeRegionMesh = ToElementMesh[WholeRegion];
WholeRegionMesh["Wireframe"]
Subscript[\[CapitalGamma], 0] =
DirichletCondition[T[r, z] == 0, {z == -Lmax, r == Rmax}];
Subscript[\[CapitalGamma], 1] =
NeumannValue[-1, {z == Lmax && r <= 1 }];
op = D[T[r, z], {r, 2}] + (1/r)*D[T[r, z], {r, 1}] +
D[T[r, z], {z, 2}];
solution =
NDSolveValue[{op == Subscript[\[CapitalGamma], 1],
Subscript[\[CapitalGamma], 0]},
T, {r, z} \[Element] WholeRegionMesh];

• Look correct to me - what would you expect as a result? – user21 Oct 14 '14 at 13:16
• Yes, the code is correct, however, it lacks some functionality, which I do not know how to apply... I would like to have dT(r,z)/dz|z=0 = alpha*[T(r,z->0+)-T(r,z->0-)] – Jerry Oct 18 '14 at 11:18

I am not exactly sure what the boundary condition is that you would like to apply but here is a way to apply a NeumannValue inside the domain (I choose a nonsense value of -1 but you get the idea):

Needs["NDSolveFEM"];
Rmax = 5; Lmax = 0.2;
WholeRegion =
ToBoundaryMesh[
"Coordinates" -> {{0, Lmax}, {0, 0}, {0, -Lmax}, {Rmax,
Lmax}, {Rmax, 0}, {Rmax, -Lmax}},
"BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {4, 5}, {5,
6}, {1, 4}, {2, 5}, {3, 6}}]}];
WholeRegionMesh = ToElementMesh[WholeRegion];
(*WholeRegionMesh["Wireframe"]*)

Subscript[\[CapitalGamma], 0] =
DirichletCondition[T[r, z] == 0, {z == -Lmax, r == Rmax}];
Subscript[\[CapitalGamma], 1] =
NeumannValue[-1, {z == Lmax && r <= 1}];
op = D[T[r, z], {r, 2}] + (1/r)*D[T[r, z], {r, 1}] +
D[T[r, z], {z, 2}];
solution =
NDSolveValue[{op ==
Subscript[\[CapitalGamma], 1] + NeumannValue[-1, z == 0],
Subscript[\[CapitalGamma], 0]},
T, {r, z} \[Element] WholeRegionMesh];

Show[
ElementMeshContourPlot[solution, PlotRange -> All,
AspectRatio -> Automatic],
WholeRegion["Wireframe"]
] 