# 3D moving heat source problem with crank-Nicholson ODE solver algorithm **URL:** \ **Category:** General\ **Created:** [June 15, 2024, 9:56pm UTC](https://fenicsproject.discourse.group/t/3d-moving-heat-source-problem-with-crank-nicholson-ode-solver-algorithm/15001 "2024-06-15T21:56:50Z")\ **Posts on this page:** 5\ **Page:** 1
**Author:** ![suparnob](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/suparnob/32/7130_2.png) [@suparnob](https://fenicsproject.discourse.group/u/suparnob)\ **Post date:** [June 15, 2024, 9:56pm UTC](https://fenicsproject.discourse.group/t/3d-moving-heat-source-problem-with-crank-nicholson-ode-solver-algorithm/15001/1 "2024-06-15T21:56:50Z")
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: ```auto 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](https://fenicsproject.discourse.group/t/newbie-trying-to-solve-a-3d-moving-heat-source-problem-with-crank-nicholson-for-a-month/14990) ---
**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\ **Post date:** [June 16, 2024, 5:08am UTC](https://fenicsproject.discourse.group/t/3d-moving-heat-source-problem-with-crank-nicholson-ode-solver-algorithm/15001/2 "2024-06-16T05:08:24Z")
> [@suparnob](#): > > ```auto > if condition.type == "Dirichlet": > bcs.append(condition.bc) > > ``` 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. ---
**Author:** ![suparnob](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/suparnob/32/7130_2.png) [@suparnob](https://fenicsproject.discourse.group/u/suparnob)\ **Post date:** [June 16, 2024, 10:51pm UTC](https://fenicsproject.discourse.group/t/3d-moving-heat-source-problem-with-crank-nicholson-ode-solver-algorithm/15001/3 "2024-06-16T22:51:11Z")
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: ```auto 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]) ---
**Author:** ![ampdes](https://avatars.discourse-cdn.com/v4/letter/a/6de8d8/32.png) [@ampdes](https://fenicsproject.discourse.group/u/ampdes)\ **Post date:** [June 21, 2024, 3:17am UTC](https://fenicsproject.discourse.group/t/3d-moving-heat-source-problem-with-crank-nicholson-ode-solver-algorithm/15001/4 "2024-06-21T03:17:42Z")
```python 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? ---
**Author:** ![suparnob](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/suparnob/32/7130_2.png) [@suparnob](https://fenicsproject.discourse.group/u/suparnob)\ **Post date:** [June 21, 2024, 3:37am UTC](https://fenicsproject.discourse.group/t/3d-moving-heat-source-problem-with-crank-nicholson-ode-solver-algorithm/15001/5 "2024-06-21T03:37:12Z")
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?