Heat flux for a moving energy source

4 months ago by
Hi all

I'm a math person and currently working on some optimization algorithms for electron beam melting. The reason I'm saying that because my backgrounds in thermo mechanics is little.

Anyway, with some help I was able to construct my simulation model on FEniCS and it runs just fine. Now I'm trying to get some results out of my model and feed them to my optimization algorithms.

Currently I need to calculate the heat flux and use them in finding the thermal gradients. I did my research and now I understand that the heat flux is the energy per unit area. Many references state that the heat flux is my governing equation (flux in = flux out). The first question is that true? and if that was the case then how can I make FEniCS report those values (excerpt of my code below)

IF not true, then how can I go about calculating that?

F = (rho*c/Dt*(T-T0)*del_T \
        + kappa*T.dx(i)*del_T.dx(i) \
        - rho*f*del_T)*dv \
        + h*(T-Ta)*del_T*da

F += sum(integrals_N)
T_M = Function(Space)
F += emis*st_bol*(T_M**4-Ta**4)*del_T*ds(5)


A = assemble(left)
b = None
T = Function(Space)

for t in numpy.arange(0, t_end, Dt):
    f.t = t
    line_n = int(abs(t / line_time))
    if (line_n % 2) == 0 and (line_n %2 < line_time):
        f.xx = (0.001 + vel*t - (length*line_n - mis))
        f.yy = 0.001 + line_n * hatch
        f.xx = (0.009 - vel*t + (length*line_n - mis))
        f.yy = 0.001 + line_n * hatch

    b = assemble(right, tensor=b)
    [bc.apply(A, b) for bc in bcs]
    solve(A, T.vector(), b, 'cg')

    timestep += 1

In case the flux was one side of my equation, I went and tried to print that value but I get the form of the function and the value.

Community: FEniCS Project

1 Answer

4 months ago by
Your form is based on the balance of internal energy without mechanical power:

$\rho\frac{\mathrm{d}u}{\mathrm{d}t}=-\nabla\cdot\mathbf{q}+\rho r$ρdudt =·q+ρr

The non-convective flux term  $\mathbf{q}$q  is called the heat flux for which you used Fourier's law

$\mathbf{q}=-\kappa\nabla T$q=κT

Now the solution field of your problem is the Temperature. So after calling solve you can use the calculated field  $T$T  to calculate the heat flux:
# Create vector valued FunctionSpace for the temperature gradient with smaller degree
VE = VectorElement('DG',Space.mesh().ufl_cell(), Space.ufl_element().degree()-1)
Vspace = FunctionSpace(Space.mesh(), VE)

# Project gradient onto function space
q = project(-kappa*grad(T), Vspace)

You can then plot that heat flux q , print it to file, print its values etc.

As for your question about the governing equation. What is that eequation supposed to govern?

Thank you Klunkean for the response

I use Gaussian function for interaction (the term 'f' in my equation). The general form (what I call governing equation) is 

$\text{ρ Cp ∂T/∂t=∇(k ∇T)+S+q}$ρ Cp ∂T/∂t=∇(k ∇T)+S+q

where S is Gaussian function, q is radiation heat loss, and there is no heat convection.
Now does that change the way you calculate the heat flux in your answer?
Thank you again for your generous response
written 4 months ago by Ashraf El Gaaly  
It doesn't. My comment on the interaction was incorrect. I only thought of the mechanical interaction and generalized, thinking of some other stresses that might act on the strains.

Of course you can model interaction through the source term (in my eq.  $\rho r$ρr ) which does not change anything regarding the heat flux.
written 4 months ago by klunkean  
I see! thank you so much
One last question and I hope that is not too much.
Back to my original question, if the heat flux is the energy per unit area, then is q (in your answer) the energy per element area for all elements?
I tried to use type(q) to find how exactly the variable look like but i didn't get any output for that command. The reason I'm asking because I want to know how to take the data and represent it.  

I really appreciate your help so much. This whole community has helped me big time, specially for someone like me who lacks the necessary background.
written 4 months ago by Ashraf El Gaaly  
It's actually power per unit area, or better put: energy per unit area per second.
In your discrete setting - assuming you use linear shape functions - you have one temperature value per node and inbetween nodes, i.e. inside the cell, you have a linear interpolation function. The gradient of temperature, and thus also the heat flux, will be a constant vector associated with a cell.
You can interpret its absolute value as the net non-convective influx (outflux if negative) of internal energy through the facets into your cell per time step.
As for visualization you could use a quiver plot with an arrow assigned to the center of each cell, or interpolate to a nodal function space.
written 4 months ago by klunkean  
I feel bad for asking so many question. please forgive me.
As you said this returns the heat flux at every node, but should there be 3 values at every node for each direction (X,Y,Z)?

If that's true, then how to return their values separately?

Excuse me if my questions sounds naive because I'm pretty naive on this 
written 4 months ago by Ashraf El Gaaly  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »