Solve a one dimensional heat transfer problem with NDSolve

Question

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.

Answer 1

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

Answer 2

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["NDSolve`FEM`"];
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 = NDSolve`VariableData[{"DependentVariables" -> {u}, 
    "Space" -> {x}, "Time" -> t}];
nr = ToNumericalRegion[mesh];
sd = NDSolve`SolutionData[{"Space" -> nr, "Time" -> 0.}];
initCoeffs = 
  InitializePDECoefficients[vd, sd, 
   "DiffusionCoefficients" -> {{-If[ElementMarker == 0, k1, k2]*
       IdentityMatrix[1]}}, "LoadCoefficients" -> {{0.}}, 
   "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.