3
$\begingroup$

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

enter image description here

where,

enter image description here

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Tim Laska Nov 6 '19 at 20:42
  • $\begingroup$ I think a time varying Neumann condition is more appropriate here to convert the point source into a flux. $\endgroup$ – Tim Laska Nov 6 '19 at 20:53
  • $\begingroup$ Version 12.1 has a tutorial on Heat Transfer Modeling $\endgroup$ – user21 Mar 20 at 19:50
6
$\begingroup$

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

Flux heating

|improve this answer|||||
$\endgroup$
  • $\begingroup$ it seems it should be $0.5g(t)$ in Neumann boundary conditions to satisfy integral heat balance. @TimLaska $\endgroup$ – Oleksii Semenov Nov 18 '19 at 10:30
  • $\begingroup$ @TimLaska.Interesting how to solve this problem in the case of absence of symmetry (different pde coefficients for example) $\endgroup$ – Oleksii Semenov Nov 18 '19 at 10:42
  • $\begingroup$ @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. $\endgroup$ – Tim Laska Nov 18 '19 at 12:58
4
$\begingroup$

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.

enter image description here

|improve this answer|||||
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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