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:

    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

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?
Is this your problem?

You are correct. That is quite a silly mistake on my end. Thanks for pointing out.

However, the main issue is the temperature. An abaqus simulation gives a solution that is one order of magnitude lower. Is my FEniCS code implementation correct?