# Correct form of source term for 3D heat equation w/ radiation

Hi everyone,

In my simulation, I wish to have a source term localized on the top surface of a cylinder. I am currently following the example shown by Bilen Emek Abali in their book ‘Computational Reality’ (pg 124). The equation of my source term is the intensity of a Gaussian laser defined as:

I(x,y,z) = \frac{2AP_{max}}{\pi w^2m}\exp{-\frac{r^2}{w^2}}\exp{-\frac{z^2}{w_2^2}}

Where A is the absorptivity of the rod, P_{Max} is the peak power of the laser, m is the mass of the cylinder, w is the width of the Gaussian beam and w_2 is the width of said beam along the z direction (I have set this to be very small such that the beam is incident on the top surface, but I am unsure if this is correct). Replicating this in Fenics, I obtain:

Pmax = 400 #peak power of laser in W
w = 0.1     #mm
R = 1.5     #mm
area = pi*R*R  #area of target (mm^2)
depth = 8.0 #thickness of target (mm)
volume = area * depth #define a disk through which the laser is absorbed
Iden = Pmax / (pi*w*w*volume*rho)    #Intensity per unit mass
absorb_depth = 1e-7

Laser = Expression('2*A*Iden*exp(-pow((x[0] - 0), 2)/(w*w)-pow((x[1]-0), 2)/(w*w)-pow((x[2]-0), 2)/(w2*w2))',degree=3, A=1, Iden=Iden,w= 0.5, w2= absorb_depth) #power (w2 localises the z-coordinates)


Abali claims in their book that the source term must be per unit mass, which I have replicated here, which seems to differ from other resources I have seen on the subject. Does anyone have any suggestions on whether or not this form is correct for a source term?

For reference, my complete code is given below:

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

parameters["allow_extrapolation"] = True
parameters["form_compiler"]["cpp_optimize"] = True
set_log_level(20)

# SOLVE THE HEAT TRANSFER PROBLEM FOR A 3D TARGET

#----------------
t_start = 0.0
t_end = 200
nstep = 200
dtime = (t_end - t_start)/ nstep  #time step
#-----------------
#Configure for dimensions of system
#----------THERMAL--PROPERTIES--[Element: W]--------
kappa = 0.0172         #conductivity [W/mm K]
c = 0.132         #Specific Heat Capacity [J/gK]
rho = 0.01925         #Density [g/mm^3]
const = kappa /(c * rho)
tau_T = 0.0 #Temperature Gradient Lag
tau_q = 0.0 #Heat Flux Lag
#---------------------------------------

pi = 3.141592653589793
T_am = 298 #ambient vacuum temperature (K)
T_a4 = T_am**4
epsilon = 0.02  # material emissivity
sigma = 5.67E-14 # W/(mm**2.K**4)
es = epsilon*sigma
Pmax = 400 #peak power of laser in W
w = 0.1     #mm
R = 1.5     #mm
area = pi*R*R  #area of target (mm^2)
depth = 8.0 #thickness of target (mm)
volume = area * depth #define a disk through which the laser is absorbed
Iden = Pmax / (pi*w*w*volume*rho)    #Intensity per unit mass
absorb_depth = 1e-7

Laser = Expression('2*A*Iden*exp(-pow((x[0] - 0), 2)/(w*w)-pow((x[1]-0), 2)/(w*w)-pow((x[2]-0), 2)/(w2*w2))',degree=3, A=1, Iden=Iden,w= 0.5, w2= absorb_depth) #power (w2 localises the z-coordinates)
#A =! emissivity for system in thermal equillibrium

geometry = mshr.Cylinder(Point(0,0,0),Point(0,0,-8),R,R) #8mm long target

# Create mesh
mesh = generate_mesh(geometry, 40)          # generate a mesh from the given geometry
Space = FunctionSpace((mesh), 'P', 1)      #define finite element function space, defined via basis functions
VectorSpace = VectorFunctionSpace(mesh, 'P', 1)
cells = MeshFunction('size_t',mesh,mesh.topology().dim()) #codimension of 0
facets = MeshFunction('size_t',mesh,mesh.topology().dim()-1) #codimension of 1
da = Measure('ds', domain=mesh, subdomain_data = facets)  #area element
dv = Measure('dx', domain=mesh, subdomain_data = cells)   #volume element

initial_T = Expression("Tini", Tini=T_am, degree=3) # extrapolate an expression for the temperature before heating
T0 = interpolate(initial_T, Space)
T = Function(Space)         # Temperature
V = TestFunction(Space)     # Test Function used for FEA
dT = TrialFunction(Space)   # Temperature Derivative
q0 = Function(VectorSpace)  # heat flux at previous time step
i = Index()
G = as_tensor(T.dx(i), (i))  #gradient of T
G0 = as_tensor(T0.dx(i), (i)) # gradient of T at previous time step

q = as_tensor(dtime/(dtime + tau_q) * (tau_q/dtime*q0[i] - kappa*(1+tau_T/dtime)*G[i] + kappa*tau_T/dtime*G0[i]),(i)) #heat
F = (rho*c/dtime*(T-T0)*V - q[i]*V.dx(i) - rho*Laser*V ) * dv + es*(T**4 - T_a4)*V*da   #final form to solve
Gain = derivative(F, T, dT)    # Gain will be usedas the Jacobian required to determine the evolution of a dynamic system

file_T = File('target3D/solution.pvd')
for t in np.arange(t_start,t_end,dtime):
print( "Time", t)
solve(F==0, T, [], J = Gain, solver_parameters={"newton_solver":{"linear_solver": "mumps", "relative_tolerance": 1e-3} }, form_compiler_parameters={"cpp_optimize": True, "representation": "quadrature","quadrature_degree": 2} )
file_T << (T,t)
q_tmp = project(q, VectorSpace)
q0.assign(q_tmp)
T0.assign(T)      #change so that T0 is equal to T for the next time step


Hi CoE2845, sorry for the mediocre answer, but unless someone else in the community is familiarized with the argument you expose, this will remain unanswered. We are surely available to solve coding questions, but seems not to be the case, but instead an issue regarding your physics. The only advice I can give you is to perform dimensional analysis, so write your equations with all the units in the variables and see how to make it match, this should solve your problem.

Best regards,
Nico

1 Like

Dear CoE2845

Sorry, I did not see this post, which is indirectly related to me, and it has been a while. Your interpretation is correct, the supply term is always a volumetric term. This is often called radiant heat (since it is mediated in a radius from a point source) and used r for its abbreviation. The typical use is a radiating heat, which is of electromagnetic origin and it goes through the material. But it is absorbed within the material and emitted only on the surface, so it is measured on the surface. This fact was the interpretation by Stefan–Boltzmann and they used T to the power of 4 on the surface.

Now, if your question is about modeling a microwave oven cooking a dish, then it is definitely this term to be used (we know that the dish is warm also inside, not only on the surface). But if you want to have a laser pointer on a surface, then it is still a volumetric term but the thickness of the source may be so small that you may model it as a surface problem. Then you can use a thickness t and multiply the power (volume) density by this t in order to find out a power surface density to be used as a boundary condition. It depends on your model. Both of these approaches are possible in FEniCS (and works, too).

Best, Emek