Hi,
I am trying to solve a PDE to simulate laser-heating of a material. The heat-source (Gaussian beam of diameter d) moves with speed v, falling on top of a thin cuboidal strip, and some heat flows out because of convection.
I am able to formulate and solve the PDE, but the values of Temperature are coming ( ~12000 K) far away from my expected experimental value (~ 1800K). Can you please look into my implementation and point any errors?
Some of my anticipations:
a) Is my solver not converging to the correct solution?
b) Am I messing up the Robin/Neumann BCs somewhere?
c) Perhaps an incorrect Variational Formulation.
d) Incorrectly setting up u_old. Do I need to use project() or interpolate() as well (apologies - I don’t udnerstand when to use these)
e) On checking min(u.vector()), I obtain -64K: which is infeasible. This could point to the possible bug.
MWE:
# Define Geometry:
length_s = 0.050
width_s = 0.020
height_s = 0.010
object_dim = np.array([length_s, width_s, height_s])
# Set physical parameters
thermal_conductivity = 7.7
T_amb = 300 # temperature of ambience
convec_coeff = 50
T0 = Constant(300.)
k = Constant(7.7) # thermal conductivity of Ti6Al4v
rho = Constant(4380) # density
cV = Constant(570)*rho # specific heat per unit volume
speed = 0.0066
laser_power = 1700
laser_diameter = 0.005
init_x = -0.40*object_dim[0]/2
init_y = 0
dt = 0.01
T = (-init_x-init_x)/speed # final time
num_steps = math.ceil(T/dt) # number of time steps
n_probes = 5
# Create Geometry
box = Box(Point(-object_dim[0]/2, -object_dim[1]/2, -object_dim[2]/2), Point(object_dim[0]/2, object_dim[1]/2, object_dim[2]/2))
mesh = generate_mesh(box, 32)
# Define expressions
laser_flux = Expression('pow(x[0]-curr_x, 2) + pow(x[1]-curr_y, 2) < pow(diameter/2, 2) ? (-8*P/(pi*pow(diameter,2)*k))*exp(-8*(pow(x[0]-curr_x, 2) + pow(x[1]-curr_y, 2))/pow(diameter,2)) : 0',
degree=2,
diameter= laser_diameter,
curr_x = init_x,
curr_y = init_y,
P = laser_power,
k = k)
r = Expression('h*k',
degree = 2,
h = convec_coeff,
k = k)
s = Expression('s',
degree = 2,
s = T_amb)
T_amb_exp = Expression('T_amb',
degree = 2,
T_amb = T_amb)
f = Constant(0) #internal heat source = 0
V = FunctionSpace(mesh, 'P', 1)
# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
u_old = Function(V)
u_old = interpolate(T_amb_exp, V)
boundary_conditions = {0: {'Neumann': laser_flux},
1: {'Robin': (r, s)},
2: {'Robin': (r, s)}}
class BoundaryTop(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[2], object_dim[2]/2, tol)
class BoundaryBot(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[2], -object_dim[2]/2, tol)
class BoundarySide(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (near(x[0], -object_dim[0]/2, tol) or
near(x[0], object_dim[0]/2, tol) or
near(x[1], -object_dim[1]/2, tol) or
near(x[1], object_dim[1]/2, tol))
# Mark boundaries
boundary_markers = MeshFunction("size_t", mesh, 2)
boundary_markers.set_all(9999)
bz1 = BoundaryTop()
bz2 = BoundaryBot()
bz3 = BoundarySide()
bz1.mark(boundary_markers, 0)
bz2.mark(boundary_markers, 1)
bz3.mark(boundary_markers, 2)
# Redefine boundary integration measure
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
# Collect Dirichlet conditions
bcs = []
for i in boundary_conditions:
if 'Dirichlet' in boundary_conditions[i]:
bc = DirichletBC(V, boundary_conditions[i]['Dirichlet'],
boundary_markers, i)
bcs.append(bc)
# Collect Neumann integrals
integrals_N = []
for i in boundary_conditions:
if 'Neumann' in boundary_conditions[i]:
if boundary_conditions[i]['Neumann'] != 0:
g = boundary_conditions[i]['Neumann']
integrals_N.append(k*g*v*ds(i))
# Simpler Robin integrals
integrals_R = []
for i in boundary_conditions:
if 'Robin' in boundary_conditions[i]:
r, s = boundary_conditions[i]['Robin']
integrals_R.append(k*r*(u - s)*v*ds(i))
F = (cV*(u-u_old)/dt*v + \
dot(k*grad(u), grad(v)))*dx + \
sum(integrals_R) + sum(integrals_N)
a, L = lhs(F), rhs(F)
# Compute solution
u = Function(V)
T_result = np.zeros((n_probes, num_steps))
x_probe = np.linspace(init_x, -init_x, n_probes)
curr_x = init_x
curr_y = 0
t = 0
for n in range(num_steps):
# Update current time
t += dt
curr_x += speed*dt
laser_flux.curr_x = curr_x
# Compute solution
solve(a == L, u, bcs)
T_result[:, n] = [u(xi, 0, object_dim[2]/2) for xi in x_probe]
# Update previous solution
u_old.assign(u)
print('%.2f'%(100*n/num_steps), '% completed')
# Completion of simulation
# Plot the results
t = np.linspace(0, T, num_steps)
plt.figure()
plt.plot(t,T_result[0,:])
plt.plot(t,T_result[1,:])
plt.plot(t,T_result[2,:])
plt.plot(t,T_result[3,:])
plt.plot(t,T_result[4,:])
plt.xlabel("t")
plt.ylabel("Temperature in K")
plt.legend(["$x={:.4f}$".format(xi) for xi in x_probe], ncol=2)
plt.show()
min(u.vector())
# gives value = - 64