How do I add a heat source inside the material?

Question

Hi guys,

I want to add a heat source inside the mesh. I have a source of radiation and a material with a specific absorption coefficient $\alpha$.
This means, the intensity $I$ of radiation which is converted into heat decays exponentially with it's depth $d$ in the material:
$I(d)=I_0 \cdot exp(-\alpha \cdot d)$
How do I implement this in my code? Or in a more mathematical way, how does this change the weak formulation?
I guess there is a better way than adding a lot of subdomains...

Best,

MaxMeier

Answer 1

If you can find a way to convert $I(d)$ into $I(x)$, where $x$ is the spatial coordinate, then you can simply define an expression with $I$ as your right hand side. As an example: suppeose that $I(x) = exp(-||x||^2)$. then your code would have to be

f = Expression('exp(-x[0]^2 - x[1]^2)')
a = ... # define bilinar form
L =f*dx

after that, you simply call the solver routines of your choice.

Two more remarks:
* I'm not 100% sure about the Expression() syntax but can't test it because it is broken on my computer
* To the best of my knowledge converting $I(d)$ into $I(x)$ for arbitrary geometries can only be achieved by solving the so-called eikonal equation. Note the part in the introduction of the wikipedia page where it says "In the special case when $F=1$, the solution gives the signed distance from $\partial \Omega$"

Comments

Ah, thank you :)
I guess the Expression syntax is correct. Only in my version of fenics I have to use x[0]*x[0] instead of x[0]**2 or x[0]^2. But this doesn't matter. The converting should be easy though, because I only want to solve carthesian or cylindrical problems.

---

I am still confused about the units: left and right hand side of the heat transfer equation gives something like $\frac{W}{m^3}$, but the intensity is $\frac{W}{m^2}$. Do I have to derivate in space or something? Or did I forget a constant?

---

Maybe you just forgot to account for the units of the constants that appear in the equation? Also, if you use the two-dimensional equation, then heat generation is not $\frac{W}{m^3}$ anymore, but $\frac{W}{m^3}$. Something similar may happen to other constants. Those are the two mistakes that I often make when checking the units in an equation or converting constants from, say, $m$ to $cm$.

---

Ah, I got it now. I have to derivate for z, due to
$\rho(u) \cdot c(u) \frac{\partial u(r,\phi,z,t)}{\partial t} = \nabla (k(u) \nabla u(r,\phi,z,t)) + \frac{\partial I(r,\phi,z,t)}{\partial z}$
when the radiation hits the surface at z = 0.