6
$\begingroup$

enter image description here

enter image description here

I'm going to solve the Laplacian equation of the electrostatic field, which consists of two triangular regions, a rectangular region, a square, and on the intersection of the two regions of $$y=x$$, there are the first and second boundary conditions. How to set the correct Boundary condition and solve the problem

tried

Ω1 = DiscretizeRegion@Triangle[{{0, 0}, {0, 1}, {1, 1}}];
(*RegionPlot[Ω1]*)
Ω2 = DiscretizeRegion@Triangle[{{0, 0}, {1, 0}, {1, 1}}];
(*RegionPlot[Ω2]*)

nv1 = NeumannValue[0, x == 0];
nv2 = NeumannValue[0, x == 1];
nv3 = NeumannValue[0, x == y];
nv4 = NeumannValue[0, x == y];

sol1 = NDSolveValue[{D[u1[x, y], x, x] + D[u1[x, y], y, y] == 
    nv1 + nv3, 
   DirichletCondition[u1[x, y] == 10, y == 1 && 0 <= x <= 1]}, 
  u1, {x, y} ∈ Ω1]
sol2 = NDSolveValue[{D[u2[x, y], x, x] + D[u2[x, y], y, y] == 
    nv2 + nv4, 
   DirichletCondition[u2[x, y] == 0, y == 0 && 0 <= x <= 1]}, 
  u2, {x, y} ∈ Ω2]

DensityPlot[sol1[x, y], {x, y} ∈ Ω1, 
 Mesh -> None, ColorFunction -> "Rainbow", PlotRange -> All, 
 PlotLegends -> Automatic]
DensityPlot[sol2[x, y], {x, y} ∈ Ω2, 
 Mesh -> None, ColorFunction -> "Rainbow", PlotRange -> All, 
 PlotLegends -> Automatic]

but the right answer should be this enter image description here

$\endgroup$
8
$\begingroup$

Your translation for the 8th equation i.e. continuity condition is wrong. This issue has been discussed in detail here so I'd like to omit corresonding explanation and simply give the solution. Values of $\gamma_1$ and $\gamma_2$ aren't given in the question so I've picked them casually.

gamma1 = 1; gamma2 = 2;
gamma = Piecewise[{{gamma1, y > x}}, gamma2];

With[{phi = phi[x, y]},
  eq = gamma Laplacian[phi, {x, y}] == 0;
  (* Alternatively: *)
  (* eq= Inactive[Div][{{gamma, 0}, {0, gamma}}.Inactive[Grad][phi,{x,y}],{x,y}] == 0; *)
  bc = {phi == 10 /. y -> 1, phi == 0 /. y -> 0};]


sol = NDSolveValue[{eq, bc}, phi, {x, 0, 1}, {y, 0, 1}];

ContourPlot[sol[x, y], {x, 0, 1}, {y, 0, 1}]~Show~Plot[x, {x, 0, 1}]

Mathematica graphics

The quality of solution above isn't that good actually, because NDSolve won't take the internal boundary at $y=x$ into consideration when discretizing the region automatically. To improve the quality, we can:

Needs["NDSolve`FEM`"]
bmesh = ToBoundaryMesh["Coordinates" -> {{0, 0}, {1, 0}, {1, 1}, {0, 1}}, 
   "BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {3, 4}, {4, 1}}], 
     LineElement[{{3, 1}}]}];
bmesh[Wireframe]
mesh = ToElementMesh[bmesh];
mesh[Wireframe]
solbetter = NDSolveValue[{eq, bc}, phi, {x, y} ∈ mesh];

Plot[{sol[x, 1/2], solbetter[x, 1/2]}, {x, 0, 1}]

Mathematica graphics

$\endgroup$
  • $\begingroup$ great answer and thanks a lot!!! $\endgroup$ – dcydhb Jul 24 '18 at 9:00
  • $\begingroup$ Very neat. As often, the problem lies in the fact that people want (i) to stick to strong formulation and (ii) to write the diffusion constants in front. $\endgroup$ – Henrik Schumacher Jul 24 '18 at 9:12
  • $\begingroup$ hi,long time no see,can you explain some code "bmesh[Wireframe] mesh = ToElementMesh[bmesh]; mesh[Wireframe]" $\endgroup$ – dcydhb Aug 9 '18 at 9:20
  • $\begingroup$ @dcydhb Please check the document of ToElementMesh and the tutorial FEMDocumentation/tutorial/ElementMeshCreation, especially the part starting from "A boundary element mesh may have internal structure; for example, to represent two material regions". $\endgroup$ – xzczd Aug 9 '18 at 10:49
  • $\begingroup$ @xzczd ok,thanks a lot! $\endgroup$ – dcydhb Aug 10 '18 at 1:33

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.