# 3D moving heat source problem with crank-Nicholson ODE solver algorithm

As advised, I am posting an abridged version of my previsouly posted code.

context: Linear transient 3D heat conduction problem.

Domain shape: Parallelopiped.

BCs: On top surface there is a moving heat source (Neumann BC). The heat source follows a line going through the middle of the top surface.

Robin BC on all other boundary faces. Ambient temp at 293 K.

initially the 3D block is at 293 K.

Problem: Temperature field becomes too high (~10^6 K).

#### Relevant parts of the code:

``````class MovingHeatSource:

self.fr = fr  # Feed rate of the moving heat source
# print(self.q)
self.t = t

def __call__(self, x):
# Calculate the current position of the moving heat source
xc = point2[0]*0.1 + self.fr * t
yc = point2[1]/2 # Assume a fixed y-coordinate for simplicity
zc = point2[2]
# print('I am here')
# Calculate the distance from the heat source center
r = ((x[0] - xc)**2 + (x[1] - yc)**2)**0.5

# Set sigma based on the torch radius for rapid decay
in_z_plane = np.isclose(x[2], zc)

sigma = self.Torch_rad / 3  # Adjust sigma as needed for faster decay
heat_expr = self.q * np.exp(-r**2 / (2 * sigma**2)) * in_z_plane

return heat_expr

#### Define domain boundaries

boundaries = [(1, lambda x: np.isclose(x[2], point2[2])),# upper z surface where heat source is moving
(2, lambda x: np.isclose(x[2], point1[2])),
(3, lambda x: np.isclose(x[1], point2[1])),
(4, lambda x: np.isclose(x[1], point1[1])),
(5, lambda x: np.isclose(x[0], point2[0])),
(6, lambda x: np.isclose(x[0], point1[0]))]

##### Assigning BC:

class BoundaryCondition():
def __init__(self, type, marker, values,dt=None):
self._type = type
if type == "Dirichlet":
u_D = Function(V)
u_D.interpolate(values)
facets = facet_tag.find(marker)
dofs = locate_dofs_topological(V, fdim, facets)
self._bc = dirichletbc(u_D, dofs)
elif type == "Neumann":
# self._bc = dt * inner(values, v) * ds(marker) #notice dt
self._bc = -inner(values, v) * ds(marker) #notice dt
elif type == "Robin":
# self._bc = values[0] * inner(T-values[1], v)* ds(marker)
self._bc = values[0] * inner(T*theta +(1-theta)*T_n -values[1], v)* ds(marker)

else:
raise TypeError("Unknown boundary condition: {0:s}".format(type))
@property
def bc(self):
return self._bc

@property
def type(self):
return self._type

##### Weak-form

F = (T * v / dt + theta * k * dot(grad(T), grad(v))) * dx - (T_n * v / dt - (1 - theta) * k * inner(grad(T_n), grad(v))) * dx

##### Solve

for i in range(num_steps):

# num_steps
print(i)

heat_source_tn = MovingHeatSource(t)

t += dt

heat_source_tnp1 = MovingHeatSource(t)

g = Function(V)
g.interpolate(lambda x: theta*heat_source_tnp1(x) + (1-theta)*heat_source_tn(x))

boundary_conditions = [BoundaryCondition("Robin", 2, (25.0, T_infty),dt=dt),
BoundaryCondition("Robin", 3, (25.0, T_infty),dt=dt),
BoundaryCondition("Robin", 4, (25.0, T_infty),dt=dt),
BoundaryCondition("Robin", 5, (25.0, T_infty),dt=dt),
BoundaryCondition("Robin", 6, (25.0, T_infty),dt=dt),
BoundaryCondition("Neumann", 1, g, dt=dt)]

for condition in boundary_conditions:
if condition.type == "Dirichlet":
bcs.append(condition.bc)
else:
F+= condition.bc

# # Define the left-hand side (LHS) and right-hand side (RHS) of the equation
L = rhs(F)
a = lhs(F)

bilinear_form = fem.form(a)
A = assemble_matrix(bilinear_form, bcs=bcs) # we write this because we do not have a dirichlet BC.
A.assemble()
solver.setOperators(A)

linear_form = form(L)
b = create_vector(linear_form)

with b.localForm() as loc_b:
loc_b.set(0)

assemble_vector(b, fem.form(L))

# Apply Dirichlet boundary condition to the vector

apply_lifting(b, [bilinear_form], bcs=[bcs])

set_bc(b, bcs=bcs)

# Solve linear problem
solver.solve(b, uh.vector)
uh.x.scatter_forward()

# Update solution at previous time step (u_n)
T_n.x.array[:] = uh.x.array

# Store the solution and the current time step
solutions.append(np.copy(uh.x.array))
time_steps.append(t)
``````

Is the solution loop correct?

1 Like

One issue I see is that this keeps accumulating boundary conditions across all time steps. Say, if you are at the 5th timestep, `bcs` now contains the `condition.bc` for timestep 1, 2, 3, 4 and 5, while I would guess that you only wanted the BC at timestep 5. Simplest workaround is to pop the last element from `bcs` before moving on to the next time.

Thanks for pointing this out. However, this piece of code actually is not being used as I donâ€™t have a â€śdirichletâ€ť bc.

I modified the following part in the solution loop:

``````    if i == 0:

boundary_conditions = [BoundaryCondition("Robin", 2, (25.0, T_infty),dt=dt),
BoundaryCondition("Robin", 3, (25.0, T_infty),dt=dt),
BoundaryCondition("Robin", 4, (25.0, T_infty),dt=dt),
BoundaryCondition("Robin", 5, (25.0, T_infty),dt=dt),
BoundaryCondition("Robin", 6, (25.0, T_infty),dt=dt),
BoundaryCondition("Neumann", 1, g, dt=dt)]

for condition in boundary_conditions:
F+= condition.bc

# F2 = F

else:
boundary_conditions_t = [BoundaryCondition("Neumann", 1, g, dt=dt)]

for condition in boundary_conditions_t:
F+= condition.bc
``````

What I do not understand is why all of sudden the temperature field becomes 0.0 after 10th timestep. The temperature is also very high:

these are max Temperature values at each time step:

array([2.95000000e+02, 1.97502723e+05, 1.58485326e+05, 3.37551641e+05,
3.55548935e+05, 5.41802417e+05, 5.44374919e+05, 6.79460856e+05,
6.86033654e+05, 8.56003231e+05, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00])

``````minsol=np.zeros(21)
for i in range(10):
minsol[i]=np.min(sol[i])
``````

Why are you computing this minsol only for 10 entries?