[Getdp] Heat equation

Matt Koch mattkoch at alum.mit.edu
Sat Dec 21 01:48:14 CET 2002


Fabio:

To begin with, isn't the heat conduction equation as follows?

rho d(cp T)/dt = +k div(grad(T)) + f(T,t)

In other words, if you set u = cp T, then your version should read as
follows!

rho d(u)/dt = +k/cp div(grad(u)) + f(T,t)

That could be the reason why you are not able to reach steady state
within the time frame you expect.

Anyway, apart from that, I received very good help from the GetDP group
on the topic of instationary heat conduction and was able to write a 3D
problem based on that input. I have attached the .geo file and the .pro
file, and it should be fairly simple to bring these from 3D to 2D.

Regards,

Matt Koch
-------------- next part --------------
// === Matt Koch ==================================== Fri-10-31-2002 ===
// Example from GetDP Mailing List!
// =====================================================================
//
// --- Group -----------------------------------------------------------
//
Group {
  SfcInr = Region[54];
  SfcOtr = Region[55];

  SfcLwr = Region[56];
  SfcUpr = Region[57];

  Vol    = Region[58];

  Sfc    = Region[{SfcInr,SfcOtr,SfcLwr,SfcUpr}];
  Dmn    = Region[{Sfc,Vol}];
}
//
// --- Function --------------------------------------------------------
//
Function {
  TheCon[Vol] = 50.0;

  TheMas[Vol] = 50.0 * 1000.0;

  TmeFct[] = ($Time < 50.0e-03) ? 1.0 : 0.0;

  PwrSrc[Vol] = 1.0e+06 * TmeFct[];

  HTCInr[SfcInr] =  5.0;
  HTCOtr[SfcOtr] = 15.0;

  TmpInr[SfcInr] = 20.0+273.15;
  TmpOtr[SfcOtr] = 30.0+273.15;

  TmpLwr[Dmn] = 10.0+273.15;
  TmpUpr[Dmn] = 12.0+273.15;

  TmeStt = 0.0;
  TmeStp = 1.0;
  TmeDel = 0.1;

  SveFct[] = !($TimeStep % 10);
}
//
// --- Constraint ------------------------------------------------------
//
Constraint {
  {Name TmpCst;
    Case {
      {Region SfcLwr; Value TmpLwr[];}
      {Region SfcUpr; Value TmpUpr[];}
    }
  }
}
//
// --- Jacobian --------------------------------------------------------
//
Jacobian {
  {Name TmpJacVol;
    Case { 
      {Region All; Jacobian Vol;}
    }
  }
  {Name TmpJacSfc;
    Case {
      {Region All; Jacobian Sur;}
    }
  }
}
//
// --- Integration -----------------------------------------------------
//
Integration {
  {Name TmpInt;
    Case { 
      {Type Gauss;
        Case { 
          {GeoElement Point      ; NumberOfPoints 1;}
          {GeoElement Line       ; NumberOfPoints 3;}
          {GeoElement Triangle   ; NumberOfPoints 4;}
          {GeoElement Quadrangle ; NumberOfPoints 4;}
          {GeoElement Tetrahedron; NumberOfPoints 4;}
          {GeoElement Hexahedron ; NumberOfPoints 6;}
          {GeoElement Prism      ; NumberOfPoints 6;}
        }
      }
    }
  }
}
//
// --- FunctionSpace ---------------------------------------------------
//
FunctionSpace {
  {Name TmpSpc; Type Form0; 
    BasisFunction {
      {Name sn; NameOfCoef Tmpn; Function BF_Node; Support Dmn; Entity NodesOf[All];}
    }
    Constraint {
      {NameOfCoef Tmpn; EntityType NodesOf; NameOfConstraint TmpCst;}
    }
  }
}
//
// --- Formulation -----------------------------------------------------
//
Formulation {
  {Name TmpFrm; Type FemEquation; 
    Quantity { 
      {Name Tmp; Type Local; NameOfSpace TmpSpc;}
    }
    Equation {
      Galerkin {    [ +TheCon[] * Dof{d Tmp}, {d Tmp} ];
                 In Vol   ; Integration TmpInt; Jacobian TmpJacVol;}

      Galerkin { Dt [ +TheMas[] * Dof{  Tmp}, {  Tmp} ];
                 In Vol   ; Integration TmpInt; Jacobian TmpJacVol;}

      Galerkin {    [ -PwrSrc[], {Tmp} ]; 
                 In Vol   ; Integration TmpInt; Jacobian TmpJacVol;}

      Galerkin { [ +HTCInr[] * Dof{Tmp}, {Tmp} ];
                 In SfcInr; Integration TmpInt; Jacobian TmpJacSfc;}

      Galerkin { [ -HTCInr[] * TmpInr[], {Tmp} ];
                 In SfcInr; Integration TmpInt; Jacobian TmpJacSfc;}

      Galerkin { [ +HTCOtr[] * Dof{Tmp}, {Tmp} ];
                 In SfcOtr; Integration TmpInt; Jacobian TmpJacSfc;}

      Galerkin { [ -HTCOtr[] * TmpOtr[], {Tmp} ];
                 In SfcOtr; Integration TmpInt; Jacobian TmpJacSfc;}

    }
  }
}
//
// --- Resolution ------------------------------------------------------
//
Resolution {
  {Name TmpRes;
    System {
      {Name Tmp; NameOfFormulation TmpFrm;}
    }
    Operation { 
      InitSolution Tmp ; SaveSolution Tmp;
      TimeLoopTheta {
	Time0  TmeStt; TimeMax TmeStp; DTime TmeDel; Theta 1.;
	Operation { 
	  Generate Tmp ; Solve Tmp;
	  Test[SveFct[]] {SaveSolution Tmp;}
        }
      }
    }
  }
}
//
// --- PostProcessing --------------------------------------------------
//
PostProcessing {
  {Name TmpPst; NameOfFormulation TmpFrm;
    Quantity {
      {Name Tmp; Value{ Local{[ {Tmp} ]; In Vol;} }}
    }
  }
}
//
// --- PostOperation ---------------------------------------------------
//
PostOperation {
  {Name TmpOps; NameOfPostProcessing TmpPst;
    Operation {
      Print[Tmp, OnElementsOf Vol, File "Tmp_3D.pos"];
    }
  }
}
//
// =====================================================================

-------------- next part --------------
// === Matt Koch ==================================== Fri-10-31-2002 ===
// Hollow Cylinder!
// =====================================================================
//
// --- Dimensions ------------------------------------------------------
//
dRdsInr =  26.0e-03;
dRdsOtr =  30.0e-03;
dHgt    = 100.0e-03;
dMsh    =   2.0e-03;
//
// --- Points ----------------------------------------------------------
//
Point(1)  = {     0.0,      0.0, 0.0, dMsh};

Point(2)  = {+dRdsInr,      0.0, 0.0, dMsh};
Point(3)  = {     0.0, +dRdsInr, 0.0, dMsh};
Point(4)  = {-dRdsInr,      0.0, 0.0, dMsh};
Point(5)  = {     0.0, -dRdsInr, 0.0, dMsh};

Point(6)  = {+dRdsOtr,      0.0, 0.0, dMsh};
Point(7)  = {     0.0, +dRdsOtr, 0.0, dMsh};
Point(8)  = {-dRdsOtr,      0.0, 0.0, dMsh};
Point(9)  = {     0.0, -dRdsOtr, 0.0, dMsh};
//
// --- Circles ---------------------------------------------------------
//
Circle(1) = {2,1,3};
Circle(2) = {3,1,4};
Circle(3) = {4,1,5};
Circle(4) = {5,1,2};

Circle(5) = {6,1,7};
Circle(6) = {7,1,8};
Circle(7) = {8,1,9};
Circle(8) = {9,1,6};
//
// --- Geometry Surfaces -----------------------------------------------
//
LopInr = newreg;
Line Loop(LopInr) = {1,2,3,4};

LopOtr = newreg;
Line Loop(LopOtr) = {5,6,7,8};

SfcGeoBtm = newreg;
Plane Surface(SfcGeoBtm) = {LopOtr,LopInr};
//
// --- Geometry Volumes ------------------------------------------------
//
VolGeo = newreg;
Extrude Surface {SfcGeoBtm, {0.0,0.0,dHgt} } {Layers{20,VolGeo,1.0}; Recombine;};
//
// --- Physics Surfaces ------------------------------------------------
//
SfcPhyInr = newreg; // 54
Physical Surface(SfcPhyInr) = {40,44,48,52};

SfcPhyOtr = newreg; // 55
Physical Surface(SfcPhyOtr) = {24,28,32,36};

SfcPhyLwr = newreg; // 56
Physical Surface(SfcPhyLwr) = {11};

SfcPhyUpr = newreg; // 57
Physical Surface(SfcPhyUpr) = {53};
//
// --- Physics Volumes -------------------------------------------------
//
VolPhy = newreg; // 58
Physical Volume(VolPhy) = {12}; 
//
// =====================================================================