# Solve a one dimensional heat transfer problem with NDSolve

The following problem is a 1-D heat transfer conduction problem:

where,

I am trying to solve with NDSolve like this way:

g[t_] := 1 + Sin[Pi*t] + Cos[Pi*t] + Sin[2*Pi*t] + Cos[2*Pi*t]

ρ = 2.7*10^3;
Cp = 1097.22;
k = 237;

NDSolve[{ρ*Cp*D[T[x, t], t] ==
D[T[x, t], {x, 2}] + g[t]*DiracDelta[x - 0.5] +
NeumannValue[0, x == 0 || x == 1], T[x, 0] == 0}, T, {x, 0,
1}, {t, 0, 1}]


It seems something is wrong with boundary condition at x direction.

• You won't be able to spatially resolve $\delta$ with a FEM mesh, but since your conductivity is constant, you can take advantage of symmetry and solve half the model converting the source into a Dirichlet condition. – Tim Laska Nov 6 '19 at 20:42
• I think a time varying Neumann condition is more appropriate here to convert the point source into a flux. – Tim Laska Nov 6 '19 at 20:53
• Version 12.1 has a tutorial on Heat Transfer Modeling – user21 Mar 20 at 19:50

Since the thermal conductivity is constant, the heat flow and temperature distributions should be symmetric about $$x=0.5$$. You should be able to convert the source term to a time dependent Neumann condition. Also, if your time scales and length scales are valid, then heat will only diffusion a small fraction of your domain in 1 second. In that case, you can simplify your model by using a Dirichlet condition at the end of the domain. Here is an example of how to convert the source to a Neumann condition using symmetry (note I think $$k$$ was inadvertently dropped in the OP). Also note that @Oleksii Semenov pointed out that the flux should be divide by 2 since it diffuses in both directions.

opts = ( Method -> {
"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.001},
"IntegrationOrder" -> 5}}});
g[t_] = 1 + Sin[Pi*t] + Cos[Pi*t] + Sin[2*Pi*t] + Cos[2*Pi*t];
rho = 2.7*10^3;
Cp = 1097.22;
k = 237;
rhocp = rho Cp;
s = NDSolve[{rhocp*D[T[t, x], t] + D[-k D[T[t, x], x], x] ==
NeumannValue[g[t] / 2, x == 0],
DirichletCondition[T[t, x] == 0, x == 1/10], T[0, x] == 0},
T, {t, 0, 4}, {x, 0, 1/10}, opts];
Plot3D[Evaluate[T[t, x] /. s], {t, 0, 4}, {x, 0, 1/20},
PlotRange -> All]
imgs = Plot[
Evaluate[If[x > 0, T[#, x] /. s, T[#, -x] /. s]], {x, -1/10,
1/10}, PlotRange -> {0, 0.00012}] & /@ Subdivide[0, 4, 100];
ListAnimate@imgs


• it seems it should be $0.5g(t)$ in Neumann boundary conditions to satisfy integral heat balance. @TimLaska – Oleksii Semenov Nov 18 '19 at 10:30
• @TimLaska.Interesting how to solve this problem in the case of absence of symmetry (different pde coefficients for example) – Oleksii Semenov Nov 18 '19 at 10:42
• @OleksiiSemenov You are correct. It was a quick and dirty answer where I did not account for symmetry. I will edit when I have time. Regarding your second question, I would probably keep it posed as a volumetric source and define a tiny interface region with average properties and create my own mesh using ToElementMesh. Here is a 2d example where the source only exists on a region. – Tim Laska Nov 18 '19 at 12:58

We can make use of low level FEM programming tools for FE matrixes calculation but load vector let's calculate "by hands". For load vector we have: $$b_i=\int_{0}^{L}g(t)\delta(x-a)\phi_i(x)dx=\begin{cases}g(t),& i=n_a\\0,& i\neq n_a\end{cases}$$ where $$\phi_i(x)$$-shape functions, $$n_a$$ - FE mesh node number which coordinate is $$a$$.

Input parameters of the model

Needs["NDSolveFEM"];
L = 0.2;                 (*dimension of calculation region*)
a = L/2;               (*position of point heat source*)
rho = 2.7*10^3;
Cp = 1097.22;
k1 = 237;                   (*conductivity in 0=<x<=a*)
k2 = 237;                   (*conductivity in a<x<=L*)
rhocp = rho Cp;
g[t_] := 1 + Sin[\[Pi]*t] + Cos[\[Pi]*t] + Sin[2 \[Pi]*t] +
Cos[2 \[Pi]*t]



Mesh generation

h = L/100;           (*spatial step size*)
Vert = Partition[Range[0., L, h], 1];
el = Table[{i, i + 1}, {i, 1, Length[Vert] - 1}];
Markers =
Table[m = el[[i]][[1]]; n = el[[i]][[2]];
If[ (Vert[[m]][[1]] + Vert[[n]][[1]])/2 < a, 0, 1], {i, 1,
Length[el]}];
mesh = ToElementMesh["Coordinates" -> Vert,
"MeshElements" -> {LineElement[el, Markers]}];
NodeSource = Floor[a/h] + 1;


FE matrices calculation

vd = NDSolveVariableData[{"DependentVariables" -> {u},
"Space" -> {x}, "Time" -> t}];
nr = ToNumericalRegion[mesh];
sd = NDSolveSolutionData[{"Space" -> nr, "Time" -> 0.}];
initCoeffs =
InitializePDECoefficients[vd, sd,
"DiffusionCoefficients" -> {{-If[ElementMarker == 0, k1, k2]*
"DampingCoefficients" -> {{rhocp}}];
methodData = InitializePDEMethodData[vd, sd];
discretePDE = DiscretizePDE[initCoeffs, methodData, sd];
{load, stiffness, damping, mass} = discretePDE["SystemMatrices"];


Solution of ODE system using Crank-Nicolson scheme

ProcDur = 4;
NumTimeStep = 400;
tau = ProcDur/NumTimeStep;   (*temporal step size*)
theta = 0.5;
xlast = ConstantArray[0, {Length[Vert]}];
X = ConstantArray[0, {NumTimeStep + 1}];
X[[1]] = {0, xlast};

Do[
time = i*tau;
A = damping/tau + theta*stiffness;
b = (damping/tau - (1 - theta)*stiffness).xlast;
b[[NodeSource]] =
b[[NodeSource]] + theta*g[time] + (1 - theta)*g[time - tau];
x = LinearSolve[A, b, Method -> "Pardiso"];
xlast = x;
X[[i + 1]] = {time, x},

{i, 1, NumTimeStep}
]; // AbsoluteTiming

ClearAll[u]
vecu = Interpolation[X];
u[t_?NumericQ] := u[t] = ElementMeshInterpolation[{mesh}, vecu[t]]

Plot3D[u[t][x], {t, 0, ProcDur}, {x, 0, L}, PlotRange -> {0, 7*10^-5},
PlotPoints -> 100]


Let's examine integral heat balance

I1 = rhocp*NIntegrate[u[ProcDur][x], {x, 0, L}, AccuracyGoal -> 8]
I2 = Integrate[g[t], {t, 0, ProcDur}]
Abs[I1 - I2]/I1
3.99644
4
0.000890264


Heat balance is fulfilled with good accuracy. Such approach also allows to set different PDE coefficient with the help of ElementMarker when symmetry relatively $$x=a$$ fails.