8
$\begingroup$

I apologize in advance if this is a somewhat basic question but I can't seem to find the solution I need and my mathematica skills are limited. I am looking to solve a 2D heat equation (for diffusion in a semiconductor) but I want to use a non-uniform mesh to improve my simulation time. I want a dense mesh at the top and a sparse mesh in lower regions of the simulations. Essentially, my question is, how do I specify the mesh for ND solve for a transient PDE?

Thanks, in advance for any help you can provide.

Here is my code with my own mesh defined:

(*Initial Parameters*)
Needs["NDSolve`FEM`"];
N0 = 1*^19;
D0 = 1.49*^-6; (*Diff Coeff for Sn in InAs*)
Q = 1.17; (*diffusion activation energy for temp dep diff*)
kBoltz = 8.617*^-5;(*Boltzmann Constant (eV)*)
T = 800 + 273; (*Diffusion Temp*)
tmax = 120*60;
xmax = 28*^-4;
ymin = -40*^-4;
meshStep = xmax*ymin/(4*^3);
Deff = D0*Exp[-Q/(kBoltz*T)];

(*Define the source Profile Conditions*)
tab = Table[Mod[.0371*x, 1], {x, 0, 28, 1}];
dcR = Round[dc, .1];
sourceTab = {};
Do[
 dc = tab[[i]];
 dcR = Round[dc, .1];
 dutyTab = Table[If[i <= dcR*10, 1, 0], {i, 1, 10, 1}];
 AppendTo[sourceTab, dutyTab],
 {i, 1, Length[tab] - 1}
 ]
sourceTab = Flatten[sourceTab];
AppendTo[sourceTab, 0];
xlist = Table[x, {x, 0, 28*10^-4, 28*^-4/(Length[sourceTab] - 1)}];
xyList = Thread[{xlist, sourceTab}];

iFun = Interpolation[xyList, InterpolationOrder -> 1, 
   Method -> "Piecewise"];

(*Create the Mesh Region*)
bmesh1 = ToBoundaryMesh[
   "Coordinates" -> {{0, 0}, {xmax, 0}, {xmax, ymin}, {0, ymin}, {0, 
      ymin/5}, {xmax, ymin/5}}, 
   "BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {3, 4}, {4, 
        1}, {5, 6}}]}];
mesh = ToElementMesh[bmesh1, 
   "RegionMarker" -> {{{xmax/2, ymin/10}, 
      1, .0001}, {{xmax/2, ymin/2}, 2, .01}}];


(*Setup and Solve Diff EQ*)
un = NDSolveValue[{
   D[u[x, y, t], t] - (D[u[x, y, t], x, x] + D[u[x, y, t], y, y])*
      Deff == NeumannValue[0, 
     y == 0 && (0 < x <= 2*xmax/5 || 3*xmax/5 <= x < xmax)],
   DirichletCondition[u[x, y, t] == N0, y == 0 && iFun[x] >= 1],
   u[x, y, 0] == 0,
   u[0, y, t] == u[xmax, y, t],
   PeriodicBoundaryCondition[u[x, y, t], x == 0, 
    TranslationTransform[{0, xmax}]]
   }, u,
  {x, 0, xmax}, {y, ymin, 0}, {t, 0, tmax},
  Method -> {"MethodOfLines", "TemporalVariable" -> t, 
    "SpatialDiscretization" -> {"FiniteElement"}}, 
  StartingStepSize -> xmax/10000
  ]

(*Plot Result*)
    Show[
     ContourPlot[un[x, y, tmax], {x, 0, xmax}, {y, ymin, 0},
      PlotPoints -> 50,
      Contours -> Flatten[Table[{1., 5., 10.}*10^n, {n, 15, 20, 1}]],
      ColorFunction -> Function[{f},
        If[f < 10^15, ColorData["Rainbow"][.0],
         ColorData["Rainbow"][(Log[10, f] - 15)/(20 - 15)]]
        ],
      ColorFunctionScaling -> False,
      PlotLegends -> BarLegend[{Function[{f},
          If[f < 10^15,
           ColorData["Rainbow"][.5],
           ColorData["Rainbow"][(Log[10, f] - 15)/(20 - 15)]]],
         {1*^15, 1*^20}}, 
        Flatten[Table[{1., 5., 10.}*10^n, {n, 15, 19, 1}]]],

      PlotRange -> All,
      ImageSize -> Large
      ], un["ElementMesh"]["Wireframe"]]
$\endgroup$
9
$\begingroup$

First initialise the problem data:

(*Initial Parameters*)Needs["NDSolve`FEM`"];
N0 = 1*^19;
D0 = 1.49*^-6;(*Diff Coeff for Sn in InAs*)Q = 1.17;(*diffusion \
activation energy for temp dep diff*)kBoltz = 8.617*^-5;(*Boltzmann \
Constant (eV)*)T = 800 + 273;(*Diffusion Temp*)tmax = 120*60;
xmax = 28*^-4;
ymin = 40*^-4;
meshStep = xmax*ymin/(4*^3);
Deff = D0*Exp[-Q/(kBoltz*T)];

(*Define the source Profile Conditions*)
tab = Table[Mod[.0371*x, 1], {x, 0, 28, 1}];
dcR = Round[dc, .1];
sourceTab = {};
Do[dc = tab[[i]];
 dcR = Round[dc, .1];
 dutyTab = Table[If[i <= dcR*10, 1, 0], {i, 1, 10, 1}];
 AppendTo[sourceTab, dutyTab], {i, 1, Length[tab] - 1}]
sourceTab = Flatten[sourceTab];
AppendTo[sourceTab, 0];
xlist = Table[x, {x, 0, 28*10^-4, 28*^-4/(Length[sourceTab] - 1)}];
xyList = Thread[{xlist, sourceTab}];

iFun = Interpolation[xyList, InterpolationOrder -> 1, 
   Method -> "Piecewise"];

Then create the finite element mesh:

mesh = ToElementMesh[
ImplicitRegion[0 <= x && 0 <= y, {x, y}], {{0, xmax}, {0, ymin}}, 
MeshRefinementFunction -> 
Function[{vertices, area}, 
Block[{x, y}, {x, y} = Mean[vertices]; 
If[y < ymin 0.1 , area > 0.000000001, area > 1]]]];
mesh["Wireframe"]

enter image description here

Solve the problem:

un = NDSolveValue[{D[u[x, y, t], 
      t] - (D[u[x, y, t], x, x] + D[u[x, y, t], y, y])*Deff == 
    NeumannValue[0, 
     y == 0 && (0 < x <= 2*xmax/5 || 3*xmax/5 <= x < xmax)], 
   DirichletCondition[u[x, y, t] == N0, y == 0 && iFun[x] >= 1], 
   u[x, y, 0] == 0, 
   PeriodicBoundaryCondition[u[x, y, t], x == 0, 
    TranslationTransform[{0, xmax}]]}, 
  u, {x, y} \[Element] mesh, {t, 0, tmax}, 
  Method -> {"MethodOfLines", "TemporalVariable" -> t, 
    "SpatialDiscretization" -> {"FiniteElement"}}, 
  StartingStepSize -> xmax/10000]

and plot the solution:

(*Plot Result*)
Show[ContourPlot[un[x, y, tmax], {x, 0, xmax}, {y, ymin, 0}, 
  PlotPoints -> 50, 
  Contours -> Flatten[Table[{1., 5., 10.}*10^n, {n, 15, 20, 1}]], 
  ColorFunction -> 
   Function[{f}, 
    If[f < 10^15, ColorData["Rainbow"][.0], 
     ColorData["Rainbow"][(Log[10, f] - 15)/(20 - 15)]]], 
  ColorFunctionScaling -> False, 
  PlotLegends -> 
   BarLegend[{Function[{f}, 
      If[f < 10^15, ColorData["Rainbow"][.5], 
       ColorData["Rainbow"][(Log[10, f] - 15)/(20 - 15)]]], {1*^15, 
      1*^20}}, Flatten[Table[{1., 5., 10.}*10^n, {n, 15, 19, 1}]]], 
  PlotRange -> All, ImageSize -> Large], 
 un["ElementMesh"]["Wireframe"]]  

enter image description here

$\endgroup$
1
  • $\begingroup$ This was a huge help, thanks a lot! $\endgroup$ – freq_show Oct 10 '17 at 22:49

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.