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.


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)/

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.


What should be the units of $Q$?


closed as off-topic by peterh, Jon Custer, user36790, ACuriousMind, Norbert Schuch Sep 16 '16 at 7:49

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "Homework-like questions should ask about a specific physics concept and show some effort to work through the problem. We want our questions to be useful to the broader community, and to future users. See our meta site for more guidance on how to edit your question to make it better" – peterh, Jon Custer, Community, ACuriousMind
If this question can be reworded to fit the rules in the help center, please edit the question.

  • $\begingroup$ Could you specify $Q(t)$ and the Neumann boundary condition specifically used? Do you mean that there's no conduction over the entire perimeter? $\endgroup$ – Gert Jul 29 '16 at 15:36
  • $\begingroup$ @Gert, does the update help? $\endgroup$ – heather Jul 29 '16 at 15:42
  • $\begingroup$ 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?' $\endgroup$ – Gert Jul 29 '16 at 15:48
  • 1
    $\begingroup$ 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... $\endgroup$ – Floris Jul 29 '16 at 19:40
  • 1
    $\begingroup$ @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. $\endgroup$ – heather Jul 29 '16 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)$.

  • $\begingroup$ 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$? $\endgroup$ – heather Jul 29 '16 at 18:40
  • $\begingroup$ 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? $\endgroup$ – heather Jul 29 '16 at 19:43
  • $\begingroup$ $\kappa=k / (\rho c)$ is the heat diffusivity ('kappa'). REe. the $\tau$, choose value as needed. $\endgroup$ – Gert Jul 29 '16 at 21:29
  • $\begingroup$ Finally, your model does not account for any convective heat losses. Is that realistic? $\endgroup$ – Gert Jul 29 '16 at 21:33
  • $\begingroup$ 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? $\endgroup$ – heather Jul 29 '16 at 23:30

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.

  • $\begingroup$ 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? $\endgroup$ – heather Jul 29 '16 at 16:48
  • $\begingroup$ 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. $\endgroup$ – user115350 Jul 29 '16 at 16:55
  • $\begingroup$ The temperature only seems to increase when the value is negative. $\endgroup$ – heather Jul 29 '16 at 16:57
  • $\begingroup$ as a diagnosis, you can use constant value for q, run simulation and see the response. $\endgroup$ – user115350 Jul 29 '16 at 17:56
  • $\begingroup$ I'm wondering if you need to define q(x, y, t). $\endgroup$ – user115350 Jul 29 '16 at 18:04

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] == 


      (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.

  • $\begingroup$ this is for a 2-d cross section of a solenoid coil. $\endgroup$ – heather Jul 29 '16 at 20:32

Not the answer you're looking for? Browse other questions tagged or ask your own question.