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$
  • $\begingroup$ Look correct to me - what would you expect as a result? $\endgroup$ – user21 Oct 14 '14 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 Oct 18 '14 at 11:18
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, privacy policy and cookie policy

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