How to increase 'stability' of Non-Linear Solution

variational formulation
1

Hi all,

I am very new to Fenics and to Finite Element Analysis in general. I am currently trying to model the heating of a 2D rectangular block in a vacuum via a laser, accounting for conduction and radiation only. Because of the nature of the problem, I do not have any Dirichlet boundary conditions, instead I only have a boundary condition which states that the energy lost to the vacuum at the boundaries must go as T**4 (ie radiative boundaries). Whilst I can obtain a solution using NonlinearVariationalProblem in order to solve:

F = T*v*dx - const*dt*dot(grad(T), grad(v))*dx + const*Q*v*dt*ds(3) - const*sum(integrals_R) -T_n*v*dx

where:

Q = Expression('2*A*t*Pmax/(pi*w*w)*exp(-2*pow(x[0]-2, 2)/(w*w))',degree=2, A=0.1, Pmax=200,w=0.3,x0=0, y0=0, t=0) #power

And sum(integrals_R) is the radiation boundary condition is defined via:

boundary_conditions = {0: {'Robin':     (es, T_a4)},
                       1: {'Robin':     (es, T_a4)},
                       2: {'Robin':     (es, T_a4)},
                       3: {'Robin':     (es, T_a4)}}


ds = Measure('ds', domain=mesh, subdomain_data = boundary_markers)
integrals_R = []
T = TrialFunction(V)
v = TestFunction(V)
T_ = Function(V)
T_n = interpolate(T_am, V)

#Account for Radiation
for i in boundary_conditions:
    if 'Robin' in boundary_conditions[i]:
        integrals_R.append(es*(T**4 - T_a4)*v*dt*ds(i)) #t-t_am will work for now to make problem solvable

I can obtain a solution using this, however, it seems wrong from the bottom face, specifically the negative temperatures. test2D

In addition, if I change any parameter (ie w in Q), the solution does not converge. Does anyone have any idea how I can make the solution more stable and hopefully fix these issues?

My full code is shown below:

from __future__ import print_function
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
from dolfin import *
from mshr import *
import math

time = 100.0            # final time
num_steps = 100     # number of time steps
dt = time / num_steps # time step size


#build rectangular mesh to model 2D slice
nx = 32
ny = 16
x0 = 0
y0 = 0
x1 = 4
y1 = 2
pi = 3.141592653589793
T_am = Constant("10") #ambient vacuum temperature
T_a4 = T_am**4
epsilon = 0.6  # material emissivity
sigma = 5.67E-8 # W/(m2.K4)
es = epsilon*sigma
Q = Expression('2*A*t*Pmax/(pi*w*w)*exp(-2*pow(x[0]-2, 2)/(w*w))',degree=2, A=0.1, Pmax=200,w=0.3,x0=0, y0=0, t=0) #power 
#----------THERMAL--PROPERTIES----------
kappa = 94         #conductivity [W/m K] 
c = 0.444          #Specific Heat Capacity [J/gK]
rho = 7            #Density [g/cm^3]
const = kappa /(c * rho)
#---------------------------------------
 
sw = Point (x0, y0)
ne = Point (x1,y1)

mesh = RectangleMesh(sw, ne, nx, ny, diagonal= "right")
V = FunctionSpace(mesh, 'P', 2)

# Define boundary condition
# build expression for power guassian (x0 and y0 displace the beam from the centre, A is the absorptivity of the target, Pmax is the maximum power of the beam, w is the width of the beam.)


#-----BOUNDARY-CONDITIONS--------------------------


#-----Define Markers for the different parts of the boundary-------
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
tol = 1E-14
class BoundaryX0(SubDomain):
    tol = 1E-14
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], x0, tol)
    
bx0 = BoundaryX0()
bx0.mark(boundary_markers, 0)

class BoundaryX1(SubDomain):
    tol = 1E-14
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], x1, tol)
    
bx1 = BoundaryX1()
bx1.mark(boundary_markers, 1)

class BoundaryY0(SubDomain):
    tol = 1E-14
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], y0, tol)
    
by0 = BoundaryY0()
by0.mark(boundary_markers, 2)

class BoundaryY1(SubDomain):
    tol = 1E-14
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], y1, tol)
    
by1 = BoundaryY1()
by1.mark(boundary_markers, 3)


# For the implimentation of these boundary conditions to be general, we can let the user specify what kind of boundary condition that applies to each of the four boundaries. We set up a Python dictionary for this purpose.


#ENCODE BOUNDARY CONDITIONS (MAY NEED ALTERING OR REMOVAL)
boundary_conditions = {0: {'Robin':     (es, T_a4)},
                       1: {'Robin':     (es, T_a4)},
                       2: {'Robin':     (es, T_a4)},
                       3: {'Robin':     (es, T_a4)}}


ds = Measure('ds', domain=mesh, subdomain_data = boundary_markers)
integrals_R = []
T = TrialFunction(V)
v = TestFunction(V)
T_ = Function(V)
T_n = interpolate(T_am, V)

#Account for Radiation
for i in boundary_conditions:
    if 'Robin' in boundary_conditions[i]:
        integrals_R.append(es*(T**4 - T_a4)*v*dt*ds(i)) #t-t_am will work for now to make problem solvable

#----------------------------------------
vtkfile = File('test2D/solution.pvd')
# fix F
# RUNTIME ERROR
F = T*v*dx - const*dt*dot(grad(T), grad(v))*dx + const*Q*v*dt*ds(3) - const*sum(integrals_R) -T_n*v*dx
T_= Function(V)   # the most recently computed solution
F = action(F,T_)
J = derivative(F, T_,T)
t = dt
bcs = []
while t <= time:
    print ('time =', t)
    Q.t = t
    problem = NonlinearVariationalProblem(F, T_, bcs, J)
    solver = NonlinearVariationalSolver(problem)
    solver.solve()
    prm = solver.parameters
    prm['newton_solver']['absolute_tolerance'] = 1E-9
    prm['newton_solver']['relative_tolerance'] = 1E-9
    prm['newton_solver']['maximum_iterations'] = 100
    prm['newton_solver']['relaxation_parameter'] = 1.0
    t += dt
    T_n.assign(T_)
    vtkfile << (T_, t)
    plot(T_)

p = plot(T_)
p.set_cmap("inferno")
#p.set_clim(0,100)
plt.colorbar(p)
plt.savefig("test2D.png")
#-----------------------------------------------------------------
2

While this isn’t the “textbook” linear heat equation, I feel like there’s a sign error in the weak form. You’d normally expect to have a weak form like: Find T s.t. \forall v

\int_\Omega (\partial_tT)v ~+ \int_\Omega c\nabla T\cdot\nabla v = \ldots\text{ ,}

while, in your code, the + from the above would be a -. Recall the change in sign due to integration-by-parts, when moving the spatial derivative from \nabla T to v.

By plugging in the test function v=T and using the product rule on the time derivative term, you can see that the sign in my equation above means that the gradient term dissipates \Vert T\Vert^2_{L^2} over time (i.e., it improves stability), while the opposite sign in your code increases this quantity (i.e., makes the system less stable).