# Laplace equation with robin boundary conditions

I want to solve the following steady state heat transfer problem with robin boundary condition at the bottom:

The following is the code for the transient solution, but how should I change the code for the steady state solution?

Clear["Global*"]
<< NDSolveFEM
L = 100;
a = 0.1;
b = 0.6;
region = Rectangle[{-L/2, -L/2}, {L/2, L/2}];
bmesh = ToElementMesh[region];

sol = NDSolveValue[{D[u[x, y, t], x, x] + D[u[x, y, t], y, y] - D[u[x, y, t], t] ==
NeumannValue[a*u[x, y, t], y == -L/2] + NeumannValue[0., x == -L/2],
u[x, L/2, t] == b,u[L/2, y, t] == 0,u[x, y, 0] == 0.
}, u, {x, y} \[Element] bmesh, {t, 0, 10},
Method -> {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement"}}, AccuracyGoal -> \[Infinity]]

DensityPlot[sol[x, y, 10], {x, y} \[Element] bmesh, PlotPoints -> 10,
MaxRecursion -> 0, Mesh -> All, ColorFunction -> "Rainbow",
PlotRange -> All, PlotLegends -> Automatic]


## Using DSolve

V 12.1 can solve this exactly.

ClearAll[u, x, y];
pde = Laplacian[u[x, y], {x, y}] == 0;
L0 = 100;
a = 1/10;
b = 6/10;
leftSide = Derivative[1, 0][u][0, y] == 0;
rightSide = u[L0, y] == 0;
bottomSide = Derivative[0, 1][u][x, 0] == a*u[x, 0];
topSide = u[x, L0] == b;
bc = {leftSide, rightSide, bottomSide, topSide};
sol = DSolve[{pde, bc}, u[x, y], {x, y}];
sol0 = sol /. K[1] -> n


(*more terms makes it more accurate*)
sol1 = u[x, y] /. First@(sol0 /. Infinity -> 200);

Plot3D[Activate[sol1], {x, 0, L0}, {y, 0, L0}]


## Using FEM

Clear["Global*"]
<< NDSolveFEM
L = 100;
a = 0.1;
b = 0.6;
region = Rectangle[{0, 0}, {L, L}]
bmesh = ToElementMesh[region]
bmesh["Wireframe"]


But triangle elements can be more accurate

bmesh = ToElementMesh[Rectangle[{0, 0}, {L, L}], "MeshElementType" -> TriangleElement]
bmesh["Wireframe"]


pde = Laplacian[u[x, y], {x, y}];
sol = NDSolveValue[{pde == NeumannValue[a*u[x, y], y == 0] ,
u[x, L] == b, u[L, y] == 0}, u, {x, y} \[Element] bmesh,
Method -> "FiniteElement"];
Plot3D[sol[x, y], {x, 0, L}, {y, 0, L}]


DensityPlot[sol[x, y], {x, y} \[Element] bmesh, PlotPoints -> 10,
MaxRecursion -> 0, Mesh -> All, ColorFunction -> "Rainbow",
PlotRange -> All, PlotLegends -> Automatic]


• A great answer. Commented Mar 28, 2020 at 23:45
• Thanks, your answer helps a lot. Commented Mar 30, 2020 at 15:28
• I just noticed that the FEM solution at the bottom is smaller than zero which is not reasonable. I tried to refine the mesh, but it does not help. Any ideas? Commented Mar 31, 2020 at 4:38
• @Jiangming I moved NeumannValue to the right side of == in the PDE specification. Now it matches OK the analytical. Also changed the mesh to triangles, I think it is more accurate, at least that what I would expect. Commented Mar 31, 2020 at 6:42
• @Nasser, you could compare the solutions for the quad and triangle mesh with the analytical solution. That would tell you which one is more accurate the this problem, given the same number of elements or better the same number of degrees of freedom (Length[mesh["Coordinates"]]) Commented Mar 31, 2020 at 7:00