Non-linear stationary heat quation in 3D

Hi there,

I’m trying to solve the heat equation in 3D without Dirichlet-BC and two Neumann-BC’s. One represents a spatially constant heat inflow that is applied to a subdomain of the boundary. The other non-linear term represents the heat radiation.

However the solver does not converge. I’m not too familiar with non-linear PDE’s, so my question is how to look for the problem.

Here you see a minimum example:

from fenics import *

def bottom(x, on_boundary):

    return x[2]<-.99 and on_boundary

def top(x, on_boundary):

    return x[2]> .99 and on_boundary

class Top(SubDomain):

    def __init__(self):

        super().__init__()

    def inside(self,x, on_boundary): 

        return top(x, on_boundary)

# define mesh and functionspace

mesh    = BoxMesh(Point(-1,-1,-1), Point(1,1,1),10,10,10)

V       = FunctionSpace(mesh, 'CG', 1)

# mark subdomain

sub_domains = MeshFunction('size_t', mesh, mesh.topology().dim()-1)

sub_domains.set_all(0)

top_domain = Top()

top_domain.mark(sub_domains, 1)

ds = Measure('ds', domain=mesh, subdomain_data=sub_domains)

# define and solve variational formulation

u0  = Constant(293.15)

u   = Function(V)

v   = TestFunction(V)

F   =   inner(grad(u) ,grad(v)) * dx + (u**4-u0**4)*v*ds - v*ds(1)       

solve(F == 0, u,  solver_parameters={"newton_solver":{"relative_tolerance": 1e-6}})

Modify your initial guess to be closer to your wanted temperature. Use for instance:

u.vector()[:] = 200

before you solve and the newton solver will converge

1 Like

Now it works, thanks!