Appropriate boundary conditions for heat equation with source [closed]

Question

I'm attempting to solve the heat equation with a time varying source term $Q$. The boundary conditions are that the temperature over the 2-D area is equal to $T_{min}$ at $t=0$ and the Neumann boundary condition is equal to $0$ to signify an insulating perimeter.

The equation that I'm using is

$\rho c \frac{\partial T}{\partial t} =k(\frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2}) + Q(t)$

The principle issue is that regardless of the value of the source term it seems to have no effect on the result. The result stays at $T_{min}$ over the entire area.

Update:

This is my code (for specifics about the equations):

rc = (8950) (385); (*density, Kg/m3 x heat capacity, J/Kg K*)
k = 385; (*thermal conductivity, W/m K*)

minT = 300;(*T, K*)
maxT = 500;(*T, K*)
time = 20;(*seconds*)

omega[t_] :=0;
q[t_] := Sin[t Pi/time];

maxX = 10/100; (*m*)
maxY = 5/100; (*m*)

fea[cell_] := {"MethodOfLines", "TemporalVariable" -> t,
  "SpatialDiscretization" -> {"FiniteElement", 
    "MeshOptions" -> {"MaxCellMeasure" -> {"Area" -> (maxX maxY)/
          cell^2}}}}

solt = NDSolve[{
    rc D[u[x, y, t], t] - 
      k (D[u[x, y, t], x, x] + D[u[x, y, t], y, y]) + q[t] == 
     NeumannValue[omega[t], True],
    DirichletCondition[u[x, y, t] == minT, t == 0]},
   u, {x, 0, maxX}, {y, 0, maxY}, {t, 0, time},
   Method -> fea[50]]

Update 2:

In light of Gert's comments:

• The boundaries are $T(x,y,0)=T_{min}$ and $\Big(\frac{\partial T}{\partial x}\Big)_{x boundary}=0$ and $\Big(\frac{\partial T}{\partial y}\Big)_{y boundary}=0$
• The shape of it is a rectangle, as shown below:

• The perimeter is perfectly insulated, so there's no conduction over the entire perimeter
• $Q(t) = \sin(t Pi/t_{max}$)...it's meant to start at zero at the beginning of the time domain, peak in the middle, and return to zero at the end of the time domain.

Bonus:

What should be the units of $Q$?

Answer 1

I second of course 'user115350's comment about $Q$, it should really be something like $Q(t)=Q_0\sin\Big(\frac{2\pi}{\tau}t\Big)$, where $\tau$ is the period of the heating cycle.

I would also recommend a small transformation like $T \to u$, where: $u=T-T_{min}$ because it makes the initial condition homogeneous:

$$u(x,y,0)=0$$ Let's have a look at the problem in one dimension first:

The PDE can then be rewritten as: $$u_t=\kappa u_{xx}+\frac{Q(t)}{\rho c}$$ Or: $$u_t-\kappa u_{xx}=\frac{Q(t)}{\rho c}$$ Solve the homogeneous PDE first: $$u_t-\kappa u_{xx}=0$$ Ansatz: $$u(x,t)=X(x)\Gamma(t)$$ Separation: $$\frac{\Gamma'}{\Gamma}=\kappa\frac{X''}{X}=-m^2$$ $$X=A\sin\Big(\frac{m}{\sqrt{\kappa}}x\Big)+B\cos\Big(\frac{m}{\sqrt{\kappa}}x\Big)$$ $$X'(0)=0 \implies B=0$$ $$X'(L)=A\sin\Big(\frac{m}{\sqrt{\kappa}}L\Big)=0$$ $$\frac{m}{\sqrt{\kappa}}L=n\pi$$ $$m^2=\frac{n^2 \pi^2 \kappa}{L^2}$$ For: $$n=1,2,3,4,...$$ $$X(x)=\sin\Big(\frac{n\pi\sqrt{\kappa}}{L}\Big)$$ $$u(x,t)=\displaystyle \sum_{n=1}^{+\infty} \Gamma_n(t)\sin\Big(\frac{n\pi\sqrt{\kappa}}{L}x\Big)$$ Now we go back to: $$u_t-\kappa u_{xx}=\frac{Q(t)}{\rho c}$$

Inserting $u(x,t)$ and expanding $\frac{Q(t)}{\rho c}$ into a Fourier series then allows to obtain (a complicated) analytical solution (with some further doodling).

But trying this approach for the two dimensional case quickly became nightmarish because its an eigenvalue problem for both $X(x)$ and $Y(x)$.

Answer 2

A few comments:

1. unit of Q will be $\frac{W}{m^3}$

2. $Q=sin(...)$ only gives you maximum heat source of 1 $\frac{W}{m^3}$. It seems too small to make significant change. Could you try $Q=1000*sin(...)$

3. My understanding Dirichlet is boundary condition but not for initial condition. I'm not sure if it is the case in your code library. You can check.

Answer 3

A few comments:

First - you have a sign error in your code:

solt = NDSolve[{
    rc D[u[x, y, t], t] - 
      k (D[u[x, y, t], x, x] + D[u[x, y, t], y, y]) + q[t] == <<<<<<<<<<
     NeumannValue[omega[t], True],
    DirichletCondition[u[x, y, t] == minT, t == 0]},
   u, {x, 0, maxX}, {y, 0, maxY}, {t, 0, time},
   Method -> fea[50]]

it should either be

      k (D[u[x, y, t], x, x] + D[u[x, y, t], y, y]) - q[t] == 

or

      (k (D[u[x, y, t], x, x] + D[u[x, y, t], y, y]) + q[t]) == 

That's why you need a negative sign on your Q in order to get a positive change in the temperature

Second - you are doing a 2D problem but using 3D units. The heat capacity of the sheet is going to be $c\rho t$ - without the thickness you are computing the heat capacity of a thick slab of material. That's OK if you really are modeling a 3D situation with symmetry, so that the heat generated is really "heat per unit length, per unit area".

Third - as I mentioned in a comment, your initial conditions didn't really allow for a very large shift in temperature - large heat capacity, small amount of heat input. The integral of your original heat input function was orders of magnitude smaller than the heat needed to get even a small change in temperature - the solver may have terminated without any shift at all due to numerical approximation issues.

Finally - you should be able to calculate the temperature change for this situation with pencil and paper. Doing so will help you figure out how to treat the units of heat input. I would recommend picking a "thin" slab (say 1 mm thick), and computing the thermal rise for a given heat per unit area.