4
$\begingroup$

I am attempting to solve the transient heat equation for a very simple geometry and I am getting a junk solution. I have identified two distinct regions (wall and super-heated fluid layer) with their own properties (density,thermal conductivity and specific heat) and solved for the steady state solution which works fine.

Steady State Solution

Then I define an initial condition for the system using part of the steady state solution.

Transient Solution Initial Condition

Finally I define the transient heat equation/boundary conditions and solve as follows:

Code Sample

The boundary conditions are basically suposed to be a constant heat flux at z = 0 and a constant temperature on the other side. Am I using the NeumannValue correctly for a constant heat flux? The results I get are junk and it throws the following two errors:

NDSolveValue::femcscd: The PDE is convection dominated and the result may not be stable. Adding artificial diffusion may help. >>
CompiledFunction::cfexe: Could not complete external evaluation; proceeding with uncompiled evaluation. >>

Can someone please help me understand exactly what the NeumannValue function is defining the way I am using it and how to fix these errors? I want the NeumannValue to be setting the first derivative of Temperature with respect to the z axis equal to a constant value.

Here is a condensed version of the code to replicate the issue:

<< "NDSolve`FEM`"

\[Rho]s = 3980;  \[Rho]sl = 958; (* kg/m3 *)
ks = .035;  ksl = .00067; (* kW/m/K *)
cs = .75; csl = 4.22; (* kJ/kg/K *)

\[Rho] = If[0 <= z < .00025, \[Rho]s, \[Rho]sl];
k =  If[0 <= z < .00025, ks, ksl];
c =  If[0 <= z < .00025, cs, csl];

eqn1 = k*\!\(
\*SubscriptBox[\(\[PartialD]\), \(z\)]\(T1[z]\)\) + 28;
Subscript[\[CapitalGamma]1, D] = 
  DirichletCondition[T1[z] == 100, z == .00048];
BCr = NDSolveValue[{eqn1 == 0, Subscript[\[CapitalGamma]1, D]}, 
   T1, {z, 0, .00048}];

Ti[z_] := \[Piecewise] {
    {BCr[z], 0 <= z < .00025},
    {100, True}
   };

eqn2 = \[Rho]*c*\!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(T[t, z]\)\) - k*\!\(
\*SubscriptBox[\(\[PartialD]\), \(z, z\)]\(T[t, z]\)\);
Subscript[\[CapitalGamma], D] = 
  DirichletCondition[T[t, z] == 100, z == .00048];
Subscript[\[CapitalGamma], N] = NeumannValue[28/k, z == 0];
soln = NDSolveValue[{eqn2 == Subscript[\[CapitalGamma], N], 
   Subscript[\[CapitalGamma], D], T[0, z] == Ti[z]}, 
  T, {t, 0, .01}, {z, 0, .00048}]
Improve this question
$\endgroup$
5
  • 2
    $\begingroup$ Could you add the code, so one can copy can paste and does not need to retype? It should also be complete such that the issue can be reproduced. $\endgroup$
    – user21
    Commented Oct 29, 2015 at 21:14
  • $\begingroup$ @user21 - I cannot comment with the entire code because it is too long. I have uploaded it using wikisend so here is a link to the code. $\endgroup$
    – dowlguest
    Commented Oct 30, 2015 at 16:17
  • $\begingroup$ wikisend.com/download/347414/1D-run.nb]1D-run.nb *here is the link to the code. If you prefer I send it to you another way just let me know. I would really appreciate your help. $\endgroup$
    – dowlguest
    Commented Oct 30, 2015 at 16:18
  • $\begingroup$ I think you should simplify your code such you can post it here. Remove everything that is not necessary to illustrate the issue. $\endgroup$
    – user21
    Commented Nov 1, 2015 at 13:24
  • $\begingroup$ @user21 - I added the condensed version as an edit on the original post. It can be copy pasted and run directly. I'm not sure why convection is an issue or how to add "artificial diffusion." $\endgroup$
    – dowlguest
    Commented Nov 2, 2015 at 16:16

1 Answer 1

Reset to default
5
$\begingroup$

When you specify that this is a temporal problem and use FEM a spatial discretization method I get to a solution:

eqn2 = \[Rho]*c*D[T[t, z], t] - k*Laplacian[T[t, z], {z}];
GD = 
    DirichletCondition[T[t, z] == 100, z == .00048];
GN = NeumannValue[28/k, z == 0];
soln = NDSolveValue[{eqn2 == GN, GD, T[0, z] == Ti[z]}, 
  T, {t, 0, .01}, {z, 0, .00048}, 
  Method -> {"MethodOfLines", 
    "SpatialDiscretization" -> {"FiniteElement"}}]

Plot3D[soln[t, z], {t, 0, 0.1}, {z, 0, 0.00048}, PlotRange -> All]

enter image description here

Does that help?

Improve this answer
$\endgroup$
2
  • $\begingroup$ That is the same solution I was getting before. I am starting to think it is because I am trying to use cylindrical coordinates rather than Cartesian coordinates. Could you tell me how to set the coordinates of my functions to be cylindrical? I tried using the SetCooridnates[Cylindrical] function but it is not working. Could this be due to being in version 10.2? I will post a seperature questions just posting this one here so you notice it. $\endgroup$
    – dowlguest
    Commented Nov 5, 2015 at 2:40
  • $\begingroup$ Laplacian takes a third argument a CoordinateChartData. Have a look there. $\endgroup$
    – user21
    Commented Nov 5, 2015 at 7:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.