# Appropriate boundary conditions for heat equation with source [closed]

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:

• 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$?

• Could you specify $Q(t)$ and the Neumann boundary condition specifically used? Do you mean that there's no conduction over the entire perimeter?
– Gert
Jul 29, 2016 at 15:36
• @Gert, does the update help? Jul 29, 2016 at 15:42
• To some it will but I'm not familiar with matlab notation. Shame, because I love this kind of problems. So far I understand this: $T(x,y,0)=T_{min}$. Also, maybe: $\Big(\frac{\partial T}{\partial x}\Big)_{x boundary}=0$ and $\Big(\frac{\partial T}{\partial y}\Big)_{y boundary}=0$. 'AmIright?'
– Gert
Jul 29, 2016 at 15:48
• The integral of your heat input over all time is (numerically) tiny - $40/\pi\approx12.7$. Also, the equation seems to imply that the heat is equally distributed over the entire area - is that correct? And the heat capacity, according to your equation, is numerically 8950*385=$3.4\cdot 10^6$. This predicts a temperature rise of $3.7\cdot10^{-6} K$. Perhaps you are doing everything right, except units... Jul 29, 2016 at 19:40
• @Floris, I've updated the $Q$ term to $\frac{-144000}{(maxX maxY maxZ)} * Sin(\frac{\pi}{\tau} t)$ which gave a 170 Kelvin rise in temperature in 4 seconds. Jul 29, 2016 at 19:48

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)$.

• I'm still looking at it, but I want $Q$ to stay in the same domain over the period, so wouldn't it be $\frac{\pi}{\tau}$ instead of $2\pi$? Jul 29, 2016 at 18:40
• So, are you saying instead of $\rho c \frac{\partial T}{\partial t} =k(\frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2}) + Q(t)$ it should be $\frac{\partial T}{\partial t} - k(\frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2}) =\frac{Q(t)}{\rho c}$? And, in yours, is $k$ the thermal conductivity? Jul 29, 2016 at 19:43
• $\kappa=k / (\rho c)$ is the heat diffusivity ('kappa'). REe. the $\tau$, choose value as needed.
– Gert
Jul 29, 2016 at 21:29
• Finally, your model does not account for any convective heat losses. Is that realistic?
– Gert
Jul 29, 2016 at 21:33
• Yeah, I understand; this was kind of a worst-case situation. Where could I find some guidance on modifying the Neumann condition to account for radiative losses? Jul 29, 2016 at 23:30

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.

• The units answer is very helpful...of course, since my area is only 2-d I can just do $\frac{W}{m^2}$. Should it be negative? Jul 29, 2016 at 16:48
• Even in 2d, there is thickness assumed. So the unit will still be $W/m^3$. In your model, you should declare the thickness assumption. The heat source can be negative, meaning it takes heat out as coolant. Jul 29, 2016 at 16:55
• The temperature only seems to increase when the value is negative. Jul 29, 2016 at 16:57
• as a diagnosis, you can use constant value for q, run simulation and see the response. Jul 29, 2016 at 17:56
• I'm wondering if you need to define q(x, y, t). Jul 29, 2016 at 18:04

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.

• this is for a 2-d cross section of a solenoid coil. Jul 29, 2016 at 20:32