Solve a one dimensional heat transfer problem with NDSolve
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
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.