6
$\begingroup$

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

enter image description here

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

Clear["Global`*"]
<< NDSolve`FEM`
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]
$\endgroup$
8
$\begingroup$

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

Mathematica graphics

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

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

Mathematica graphics

Using FEM

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

Mathematica graphics

But triangle elements can be more accurate

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

Mathematica graphics

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}]

Mathematica graphics

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

Mathematica graphics

| improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ A great answer. $\endgroup$ – Please Correct GrammarMistakes Mar 28 at 23:45
  • $\begingroup$ Thanks, your answer helps a lot. $\endgroup$ – Jiangming Mar 30 at 15:28
  • $\begingroup$ 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? $\endgroup$ – Jiangming Mar 31 at 4:38
  • $\begingroup$ @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. $\endgroup$ – Nasser Mar 31 at 6:42
  • 2
    $\begingroup$ @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"]]) $\endgroup$ – user21 Mar 31 at 7:00

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.