2
$\begingroup$

I am trying to solve simple 3D Poisson's equation of thermal conduction due to heat generation (q) in a spherical composite, where from 0 to a (inner radius) conductivity is k1 and from a to b (outer radius) conductivity is k2. And BC is temperature becomes zero at r=b. I am using FEM to solve this problem, but the result of NDSolveValue is coming only till b=4.99 not b=5. Please help. I guess the output is also not correct.

  Needs["NDSolve`FEM`"]
    k1=10;k2=1;a=1;q=1;b=5;
    kk[x_, y_, z_] := 
      Boole[x^2 + y^2 + z^2 < a^2] k1 + Boole[x^2 + y^2 + z^2 >= a^2] k2;
    area = ImplicitRegion[x^2 + y^2 + z^2 - b^2 <= 0, {x, y, z}];
    bc1 = DirichletCondition[t[x, y, z] == 0, x^2 + y^2 + z^2 - b^2 == 0];
    op = kk[x, y, z] Laplacian[t[x, y, z], {x, y, z}] + q;
    tfun = NDSolveValue[{op == 0, bc1}, t, {x, y, z} ∈ area];
    Plot[tfun[x, 0, 0], {x, 0, 5}]
$\endgroup$
4
  • $\begingroup$ tfun[x,y,z] is evaluable only, if the point x,y,z lies inside the grid(mesh) of your FEM-solution. Obviously point 5,0,0lies outside! $\endgroup$ Jun 21, 2020 at 21:13
  • $\begingroup$ Clarification: why do you think output is not correct? Also variable "area" is misleading. Domain is a volume. Like the use of Boole. $\endgroup$ Jun 21, 2020 at 23:10
  • $\begingroup$ Thank you for replying. @Ulrich Neumann Sorry I didn't understand why 5,0,0 lies outside the mesh. $\endgroup$
    – Nilabh
    Jun 22, 2020 at 4:16
  • $\begingroup$ Thank you for replying. @PaulCommentary I think output is not correct because when I am trying to solve by DSolve by assuming temperature and heat flux continuity at the interface, it is giving me different result. $\endgroup$
    – Nilabh
    Jun 22, 2020 at 4:19

1 Answer 1

3
$\begingroup$

When a symbolic region is meshed there can be discrepancies between what the symbolic region represents and what the numerical approximation represent; that is especially the case for curved regions. "OpenCascade" typically does a good job at preserving the boundary. So if you can, it's advisable to represent the region with graphics primitives such that "OpenCascade" can work like so:

Needs["NDSolve`FEM`"]
k1 = 10; k2 = 1; a = 1; q = 1; b = 5;
kk[x_, y_, z_] := 
  Boole[x^2 + y^2 + z^2 < a^2] k1 + Boole[x^2 + y^2 + z^2 >= a^2] k2;
area = Ball[{0, 0, 0}, b];
bc1 = DirichletCondition[t[x, y, z] == 0, True];
op = kk[x, y, z] Laplacian[t[x, y, z], {x, y, z}] + q;
mesh = ToElementMesh[area, "BoundaryMeshGenerator" -> "OpenCascade"];
tfun = NDSolveValue[{op == 0, bc1}, t, {x, y, z} \[Element] mesh];
(*Plot[tfun[x,0,0],{x,0,5}]*)
tfun[b, 0, 0]

(* 1.13485*10^-16 *)

Other options would be to refine the mesh. Such that the point {b,0,0} is actually part of the mesh or change the interpolating function property such that extrapolation will be performed.

$\endgroup$
2
  • $\begingroup$ Thank you so much. $\endgroup$
    – Nilabh
    Jun 22, 2020 at 14:36
  • $\begingroup$ @Nilabh, you are welcome. $\endgroup$
    – user21
    Jun 22, 2020 at 19:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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