Converting complex TD Heat Equation to Time Independent form

variational formulation
1

Hi everyone,

I have a complex heat equation that takes into account radiation into a vacuum (T_{0} = 0) whilst the top surface of my block is being heated by a laser source with a guassian profile:
Q = \frac{2AP_{max}}{\pi w^2}\exp(-2\frac{r^2}{w^2})
Where A is the absorptance of the block, P_{max} is the peak power of the laser and w is the FWHM of the guassian beam. When I have previously solved this problem in 2D rectangular coordinates (I am now trying to solve for 3D cylindrical coordinates), I used:

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_)

However, I now wish to find the steady state solution of this problem, wherein radiative cooling and heating via the laser balance themselves out across the whole mesh.

In addition, I also wish to convert this problem into cylindrical coordinates but I am unsure how I can convert (x,y,z) =(x[0],x[1],x[2]) into the required cylindrical coordinates.

My previous 2D time-dependent rectangular code is shown below. Thank you in advance

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 = 180.0            # final time
num_steps = 180     # number of time steps
dt = time / num_steps # time step size

#Configure for dimensions of system
#build rectangular mesh to model 2D slice
nx = 48
ny = 16
x0 = 0
y0 = 0
x1 = 12  #mm
y1 = 4    #mm
pi = 3.141592653589793
T_am = Constant("0") #ambient vacuum temperature
T_a4 = T_am**4
epsilon = 0.05  # 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]-6, 2)/(w*w))',degree=2, A=0.1, Pmax=500,w=0.5, t=0) #power 
#----------THERMAL--PROPERTIES----------
kappa = 0.174         #conductivity [W/mm K] 
c = 0.132         #Specific Heat Capacity [J/gK]
rho = 0.019            #Density [g/mm^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  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')
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")
plt.savefig("test2D.pdf")