0

I'm trying to solve the 2D axisymmetric heat equation, where only one section of the tube is being heated up. I have written the following code below:

q=Piecewise[{{0,z<=0},{10000t,0<=z}}]
bc1=q==1.38(D[T[r,z,t],r]/.r->0.03);
bc2=(D[T[r,z,t],r]/.r->0.00001)==0;
ic1=T[r,z,0]==273.15;

heat=(1521740*D[T[r,z,t],t])==D[T[r,z,t],r,r]+D[T[r,z,t],r]+D[T[r,z,t],z,z];
sol=NDSolve[{heat,bc1,bc2,ic1},{T},{r,0.00001,0.03},{z,-0.2,0.2},{t,0,100}]

However, I get the error NDSolve: Encountered non-numerical value for a derivative at t==0`. The code works when I change q=10000t instead. I think something is likely wrong with the piecewise boundary condition but I'm not sure what is wrong.

Edit: After updating q, the code is able to run but the result seems unphysical (The temperature in some section decreases). It seems like the boundary condition bc1 isnt being properly obeyed.

q=Simplify`PWToUnitStep@Piecewise[{{0,z<=0},{10000t,0<=z}}]
bc1=q==1.38(D[T[r,z,t],r]/.r->0.03);
bc2=(D[T[r,z,t],r]/.r->0.00001)==0;
bc3 = (D[T[r, z, t], z] /. z -> 0.2) == 0;
bc4 = (D[T[r, z, t], z] /. z -> -0.2) == 0;
ic1=T[r,z,0]==273.15;

heat=(1521740*D[T[r,z,t],t])==D[T[r,z,t],r,r]+D[T[r,z,t],r]+D[T[r,z,t],z,z];
sol=NDSolve[{heat,bc1,bc2,bc3,bc4,ic1},{T},{r,0.00001,0.03},{z,-0.2,0.2},{t,0,100}]
Interpolation[Table[{r, sol[[1, 1, 2]][r, -0.01, 10]}, {r, 0.0001, 0.03, 
 0.0001}]]'[0.03]

The interpolation helps to evaluate the partial derivative of r at z=-0.01. Based on the Piecewise function bc1, this should be 0. However, the interpolation gives me -6247.

user21
  • 39,710
  • 8
  • 110
  • 167
Lucas
  • 103
  • 5
  • To be more specific, change definition of q to q = Simplify`PWToUnitStep@Piecewise[{{0, z <= 0}, {10000 t, 0 <= z}}] avoids the ndnum warning. You still need to fix the bcart warning, though. Please notice bcart is a serious problem: https://mathematica.stackexchange.com/q/73961/1871 – xzczd Mar 07 '23 at 09:43
  • Instead of bc1 and bc2, include a 10000 t /1.38 NeumannValue[1, r == .03 && z <= 0 ] in your pde.. – Ulrich Neumann Mar 07 '23 at 10:24
  • @xzczd Introducing NeumannValue simply solves the problem I think. No need to use the link! That's why I propose to reopen the question – Ulrich Neumann Mar 07 '23 at 10:48
  • @UlrichNeumann 1. Since OP wrote "I'm not sure what is wrong" I think the topic of the post is to understand what's wrong with the code. If OP clarifies further and make the post not a duplicate, I'll vote to reopen. 2. I think turning to FEM isn't a good idea in this case, because FEM will automatically add zero Neumann condition in z direction, but at the moment we just don't know what b.c.s OP needs in z direction. – xzczd Mar 07 '23 at 11:21
  • @xzczd Ok let's wait, thanks. But it's not to difficult to add convection bc[z] with FEM... – Ulrich Neumann Mar 07 '23 at 11:35
  • @xzczd Yes changing the definition of q helped to solve the issue. But I'm not entirely sure how it works... does q get changed into a single function? As for the bcart warning, I forgot to add the convection boundary conditions at the two ends. – Lucas Mar 08 '23 at 10:07
  • In short, your original q doesn't work because of a bug of NDSolve, it's not your fault. (This might be related to the improper differentiation of Piecewise[…]. ) – xzczd Mar 08 '23 at 12:59
  • I see thanks, I've edited the code now and it's able to run, but it seems like the BC may not be properly applied. I've edited the question to add my new code and to describe the issue that arises. Thanks for your patience :) – Lucas Mar 08 '23 at 13:38
  • Well, since the problem has totally changed, you should post a new question for it. Anyway, what you've observed is not surprising, because the Piecewise function is approximated with finite difference formula in this case. (If you choose FEM for spatial discretization, the situ will be slightly different, but similar. ) Related: https://mathematica.stackexchange.com/q/10055/1871 – xzczd Mar 09 '23 at 05:59

0 Answers0