Point Source moving

4 months ago by
Hello everyone,

I am solving the heat equation with a moving point source in a solid. I am wondering how the heat source is interpolated among nodes.

Suppose the solid box is 1.0*1.0*0.4. Point source moves along x direction. Boundary conditions: bottom Dirichlet, sides and top are adiabatic.

I did two simulations: (1) point source moving along x direction at z=0.375 plane, I set the conductivity of nodes above the scan plane (z>0.375) to be zero (2) remove the elements above scan plane and point source moving in the same manner as (1). Everything else (including mesh) is the same for this two simulations.

Here is a schematic:

Then I extracted out the temperature history (a short time) of a same point from the two simulations. I expect them to be same since zero conductivity is equal to adiabatic bcs. But the results are quite different:

Thought: from the results, it seems some heat loss to the nodes above the scan layer. So I extracted out temperature for the dofs above the scan layer, they are all zeros.

Not sure where the discrepancy comes from. Any help will be greatly appreciated!

Is it possible that the point source is being extrapolated to the nodes but not with the FE shape functions?

Below is the script for simulation 1. For simulation 2, I just changed zmax to be 0.375 and kept everything else to be the same.

from dolfin import *
import time
import numpy as np

#define domain
x_min = y_min = z_min = 0.0
x_max = y_max = 1.0
z_max = 0.4 

#define number of elements
nx = ny = 40 
nz = 16
#create mesh
mesh = BoxMesh(Point(x_min, y_min, z_min), Point(x_max, y_max, z_max), nx, ny, nz)

#define point source characteristics
velocity_array = np.array([491,0,0],'d')

# Time-stepping parameters
t_end = 2/491 
dt = (1./30.)*(hmin)**2*0.25 # factor is empirical max step size

# Initial condition and body forces 
ic = Constant(0.0)
f = Constant(0.0)

# define function space
V = FunctionSpace(mesh, 'P', 1)

# Create boundary markers
def boundary(x, on_boundary):
    return near(x[2],z_min) # Define boundary condition
bc = DirichletBC(V, Constant(0.0), boundary)

# Define trial and test function and solution at previous time-step
u = TrialFunction(V)
v = TestFunction(V)
u0 = Function(V)

# apply differenet material properties kappa and cp

# rhs (note that the sign of the f term might be wrong) 
L=-dt*(kappa*dot(grad(u0), grad(v))*dx - f*v*dx)

# lhs

# Prepare solution function and solver
u = Function(V)

# assemble mass matrix
A = assemble(a)

#lump mass matrix
unity = Function(V)
unity.vector()[:] = 1.0 
A_lumped = A*unity.vector()  # lump mass via row sum

# Prepare initial condition

dof_x=V.tabulate_dof_coordinates().reshape(V.dim(), -1)

z = dof_x[:, 2]
indices = np.where(np.logical_and(z > 0.375, z <= 0.4))[0]
# Get coordinates of dof
xs = dof_x[indices]

# Time-stepping
t = 0.0
while t < t_end:

    b = assemble(L)

    if t*velocity_array[0]<end_distance :
       delta = PointSource(V, Point(0.5+t*velocity_array[0], 0.5+t*velocity_array[1], 0.375+t*velocity_array[2]), dt)
    # Euler explicit
    u.vector()[:] = u0.vector() + b.array()/A_lumped.array()
    # printing the temperature at a point
    if step_number%5==0:
        print(t,' ',u(0.5,0.5,0.375-.06),' ','\n',file=open('moving1.txt','a')) 
    # Move to next time step
    t += dt

# extract out the dofs above the scan plane
vals = u.vector()[indices]
for z, v in zip(xs, vals):
    print(z,' ',v,file=open('moving1_vertice.txt','a'))
Community: FEniCS Project

1 Answer

4 months ago by
For nodes on the z=0.375 boundary plane in case 1, there is no zero-flux condition, and even otherwise we already expect the linear elements to decay to 0 on the adjoining nodes (z>0.375, where conductivity is zero).  Thus, there is additional heat in the volume outside your quasi-boundary: where it decays from nodes @ z=0.375 to the next nodes at z>0.375  In case 2, there are no adjoining nodes and zero flux is asserted, so the nodes have no decay. 

This is just a guess at the issue.  If you compute an accurate integral of the heat over the entire volume for each case, I would expect them to be near the same, thus the inner nodes must decrease their values to account for the heat on the quasi-boundary decay volume in case 1 vs. case 2?  Also, if you used a finer mesh, I would expect the error-difference to decrease, which you could check with a convergence test.
Thanks for the information, jwinkle!
If I set the thermal conductivity to be zero for all the nodes above z=0.375, even if there is no zero-flux condition, I do not see how heat will be lost/decay to nodes above this plane. Also I output the temperature fields for all the dofs above z=0.375, they are all zeros which kind of show that no heat loss at all to the nodes above.

"an accurate integral of the heat over the entire volume for each case"-- good idea. I'll try this!

Convergence test: i ran simulations with finer mesh already, the difference is still there.
written 4 months ago by Peipei  
Yes, you are correct that the dofs would be zero above the boundary plane.  Yet still, the elements will in case 1 be interpolated (linearly, with your degree = 1 elements) to zero to these adjoining dofs, as this is the definition of the linear basis functions per se.  I really don't know if this at all makes a difference to the solution, it is just the only difference I can see between the two cases.   That is, in case 1, the FE solution will decay from values at z=0.375 to zero at the next dofs, so there is non-zero (interpolated) heat just outside z=0.375.  In case 2, it is a hard-stop as there are no dofs at all outside that plane.

And if in fact with finer meshes, the error is numerically the same between the cases, it's probably not due to this anyway, but it would be interesting to see if there is a difference in the integral of the solutions.
written 4 months ago by jwinkle  
Thanks a lot, @jwinkle !

I think I sort of get it. Even the dofs above z=0.375 return zeros. There will be heat loss due to the linear interpolation per se.

I just did the integral of solutions using assemble() in fenics. Case 1 and case 2 returned similar values. Then for case 1, I assembled the solution over a subdomain(z<=0.375), it returned a much smaller value comparing to the previous ones. So sort of proved what you said.

Again, thanks for all the help!

PS: I am simulating AM process using the quiet FEM (i.e. set zero thermal conductivity to the elements that haven't been solidified.) But apparently this will give me certain amount of errors like what I showed above. Do you know if I can set the connectivity of nodes above the current scan layer to be zeros in fenics so it can "pretend" to be a hard-stop?
written 4 months ago by Peipei  
@Peipei:  perhaps you can use a discontinuous Galerkin method in this case, which allows step changes across the dofs, but I am by no means an expert in this.  --JW
written 4 months ago by jwinkle  
Got it. I'll take a look at how discontinuous Galerkin method works.

Thank you! :)
written 4 months ago by Peipei  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »