3
$\begingroup$

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["NDSolve`FEM`"];
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];
$\endgroup$
2
  • $\begingroup$ Look correct to me - what would you expect as a result? $\endgroup$
    – user21
    Commented Oct 14, 2014 at 13:16
  • $\begingroup$ 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-)] $\endgroup$
    – Jerry
    Commented Oct 18, 2014 at 11:18

1 Answer 1

3
$\begingroup$

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["NDSolve`FEM`"];
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"]
 ]

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.