Point Source moving
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 hatch_spacing=225.0/5000 velocity_array = np.array([491,0,0],'d') end_distance=0.5-hatch_spacing/2 # Time-stepping parameters t_end = 2/491 hmin=mesh.hmin() dt = (1./30.)*(hmin)**2*0.25 # factor is empirical max step size total_number_of_steps=t_end/dt # 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,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 string0='x<=0.375+tol?k_0:k_1' kappa=Expression(string0,degree=0,tol=1e-4,k_0=1.0,k_1=0.0) # 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 a=u*v*dx # Prepare solution function and solver u = Function(V) # assemble mass matrix A = assemble(a) bc.apply(A) #lump mass matrix unity = Function(V) unity.vector()[:] = 1.0 A_lumped = A*unity.vector() # lump mass via row sum # Prepare initial condition u0.interpolate(ic) u.interpolate(ic) 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)) # Get coordinates of dof xs = dof_x[indices] # Time-stepping t = 0.0 step_number=1 while t < t_end: b = assemble(L) if t*velocity_array<end_distance : delta = PointSource(V, Point(0.5+t*velocity_array, 0.5+t*velocity_array, 0.375+t*velocity_array), dt) delta.apply(b) bc.apply(b) # 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 u0.assign(u) t += dt step_number+=1 # 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'))
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.