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:
def __init__(self,t, Q=100.0, Torch_rad=Torch_rad, fr=fr):
self.Torch_rad = Torch_rad # Radius of the heat source
self.fr = fr # Feed rate of the moving heat source
self.q = Q/(np.pi*Torch_rad**2)
# 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)
# being an implicit we start with t+dt
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])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
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?
Previous post:
Newbie trying to solve a 3D moving heat source problem with crank-Nicholson for a month - General - FEniCS Project